Finding the traffic intensity, A in Erlang B

I am given the 1) blocking probability, B

2) Number of Trunk , N

Erlang B Equation

B = (A^N/factorial(N))/sum(A^I/factorial(I),I,0,N)

How do i code to find the Traffic intensity, RHO in erlang B equation?

%%ACell = the RHO i want to find which is the A
%%N = trunk
%%FN = Factorial trunk
%%%B = blocking probability
syms I A
B = 0.05;
N = 55;
FNoC = factorial(N);
ACell = solve((A^N)/FN == B*symsum(A^I/factorial(I),I,0,N),A);

Answer 1

[David Wilson]

OK, that's some nasty equation you've got there. But the symbolic toolbox did find some solutions, in fact it found many (& complex)

syms I A
B = 0.05;
N = 55;
FN = factorial(N);
ACell = solve((A^N)/FN == B*symsum(A^I/factorial(I),I,0,N),A);

It returned 55 of them, 54 complex and 1 real.

>> length(ACell)
ans =
    55

Let's convert to double, and take a look at the only real one:

    >> A = double(ACell);
    >> A = A(end)
    A =
  49.539390342298994

Seems to indicate that A is a little over 49.5. Does this make sense?

>> N = 55; 
>> Erlang = @(A) (A^N/factorial(N))/sum(A.^([0:N])./factorial([0:N])) % Very ugly

>> Bx = Erlang(49.539390342298994) % should return original B = 0.05

Bx =
   0.050000000000000

Of course the large factorials and ugly expression could be made far more robust, but it seems to work. For example probably should use cumprod instead of factorial, e.g:

>> factorial([1:7])
ans =
           1           2           6          24         120         720        5040
>> cumprod([1:7])
ans =
           1           2           6          24         120         720        5040

[Kenny Kwan]

Yes, It helps alot. Thank you for your help and your time .

But it is weird how when i follow your steps. Matlab shows a different value of A(end) for me.

>> syms I A
>> B = 0.05;
>> N = 55;
>> FN = factorial(N);
>> ACell = solve((A^N)/FN == B*symsum(A^I/factorial(I),I,0,N),A);
>> A = double(ACell);
>> A = A(end)
A =

  36.5271 +14.2925i

There's no real number when i type in the A = double(ACell)

[Yousry AbdelHamid]

It works absolutely perfect with me. I am using R2022b. Thank you so much. Could you please help with some detail the use of "double" format here? Thanks again.

Answer 2

[David Wilson]

OK, perhaps there is no guarentee that the symbolic returns the real solution last, so in that case you need to resort to a numerical strategy, perhaps using fsolve:

>> N = 55; B = 0.05; 
>> Erlang = @(A) (A^N/factorial(N))/sum(A.^([0:N])./cumprod([0,0:N-1]+1)) % Very ugly

>> A = fsolve(@(A) Erlang(A)-B, 40) % start guess is say A = 40
A =
   49.5409
>> Bx = Erlang(A) 
Bx =
    0.0500

The function could be improved.

[Kenny Kwan]

Alright, Thank you for your help. I'll try on it and work it out. Have a nice day my friend.

Answer 3

[Nermin Tawfik]

I want the value of A at N=197 & B=0.05