How to solve 'n' number system of equations

What I am trying to do here is find PL as equation of variables, then differentiate it to find dpl/dpi ,subtitute it in this formula

where Li is written as k in my code then multiply it with df(i) and equate it to L which is 10.

The issue is if i length of 'a' is 2 then this code will work properly but if length of 'a' is greater than 2, this will cause for k(i)*df(i) to have more than one variable inside ( p1, p2 ,p3 ...) and i cant solve the system. I tried to convert it to matrix and solve it but it was not possible due to it being non linear.

clc
clear
a=[0.012;0.018;0.014];%The cofficeient next to power in the IFC equation 
ng =length(a);
syms p [1 ng] 
symvar Y
L=10;%Value of lambda
b=[7,6.3,6.8];%The cofficeientin the IFC equation
B=[0.1 0.01 0.01; 0.01 0.13 0.02; 0.01 0.02 0.04];
for i=1:length(a)
    for j=1:length(a)
        pl(j,i) = p(i)*B(i,j)*p(j);
    end
    Pl(i) = sum(pl(:,i));
end
Pl = sum(Pl);
for i=1:length(a)
df(i)=a(i)*p(i)+b(i);
dp(i)=diff(Pl,p(i));
k(i)=1./(1-dp(i));
      Y(i)=linsolve(k(i)*df(i)==L,p(i))
      plc=subs(pl(i),Y(i));
n(i)=sum(plc);
end
pt=sum(Y)-plc;
disp('The value of Power of each plant respectively in a matrix are')
disp((Y))
fprintf('The value of power losses is = %f ',abs(plc))
fprintf('\nThe value of total power is = %f  ', abs(pt))

Answer 1

[Matt J]

It appears to me that the equation can be re-arranged quite easily into a set of linear equations in P,

or

which is easily solved in one line:

a=[0.012;0.018;0.014];%The cofficeient next to power in the IFC equation 
b=[7,6.3,6.8];
L=10;%Value of lambda

B=[0.1 0.01 0.01; 0.01 0.13 0.02; 0.01 0.02 0.04];

P=(B+diag(a)/(2*L))  \ (1 - b(:)/L   )/2

P = 3×1

1.0847 0.8331 3.2553

[Osama Aliwat]

I am not sure about that.

To put it in simple terms, after multiplying k with df I will end up with three equations:

x1P1+y1P2+z1P3 EQ1

x2P1+y2P2+z2P3 EQ2

x3P1+y3P2+z3P3 EQ3

where x,y,z is some random number. and I equate them to L which is 10 and i will have 3 equations with 3 unknowns which I want to solve.

x1P1+y1P2+z1P3=L

x2P1+y2P2+z2P3=L

x3P1+y3P2+z3P3=L

and I am getting these equations as a matrix Y 1x3

[Matt J]

Well, if you don't believe me, you can go ahead an substitute the values I'm gettin for P into your equations and see if the equations are satisfied.

[Osama Aliwat]

I am sorry but thats just not correct. Even after hand calculations the answer is not correct.

We can find the three equations from k(i)*df(i) and equate it to L

using this code again

clc
clear
% a=[0.01;0.015]
a=[0.012;0.018;0.014];%The cofficeient next to power in the IFC equation 
ng =length(a);
syms p [1 ng] 
L=10;%Value of lambda
% b=[8.5;9.5]
% B=[4e-4 0; 0 0]
b=[7,6.3,6.8];%The cofficeientin the IFC equation
B=[0.1 0.01 0.01; 0.01 0.13 0.02; 0.01 0.02 0.04];
for i=1:length(a)
    for j=1:length(a)
        pl(j,i) = p(i)*B(i,j)*p(j);
    end
    Pl(i) = sum(pl(:,i));
end
Pl = sum(Pl);
for i=1:length(a)
df(i)=a(i)*p(i)+b(i);
dp(i)=diff(Pl,p(i));
k(i)=1./(1-dp(i));
       Y(i)=linsolve(k(i)*df(i)==L,p(i));
      X(i)=simplify(k(i)*df(i)==L,'ignoreanalyticconstraints',true)
      plc=subs(pl(i),Y(i));
n(i)=sum(plc);
end
[A,B]=equationsToMatrix(X,p );
plc=sum(n)
pt=sum(Y)-plc;
disp('The value of Power of each plant respectively in a matrix are')
disp((Y))
fprintf('The value of power losses is = %f ',abs(plc))
fprintf('\nThe value of total power is = %f  ', abs(pt))

X find the equations but i cannot take the equationstomatrix because it shows me vars are not linear and if we run the code we can see that X gives us

X =
 
[503*p1 + 50*p2 + 50*p3 == 750 & 3*p2 + 3*p3 ~= 17650,
    100*p1 + 1309*p2 + 200*p3 == 1850 & p2 ~= -350,
    100*p1 + 200*p2 + 407*p3 == 1600 & p3 ~= -3400/7]

503*p1 + 50*p2 + 50*p3 == 750

100*p1 + 1309*p2 + 200*p3 == 1850

100*p1 + 200*p2 + 407*p3 == 1600

which if we solve, we will not get the same result as said above.

Now I am trying to get to solve these equations from X which I could not figure it out.

[Matt J]

Isn't the equation you ultimately want to satisfy k.*df=10? The code below demonstrates that this is achieved:

a=[0.012;0.018;0.014];%The cofficeient next to power in the IFC equation 
b=[7;6.3;6.8];
L=10;%Value of lambda

B=[0.1 0.01 0.01; 0.01 0.13 0.02; 0.01 0.02 0.04];

P=(2*B + diag(a/L))  \ (1 - b./L )

P = 3×1

1.0847 0.8331 3.2553

dp=2*B*P;
df=a.*P+b;
k=1./(1-dp);

k.*df  %Check the original equation

ans = 3×1

10.0000 10.0000 10.0000

We can also verify the equations you presented, though they still are somewhat obscure to me:

[p1,p2,p3]=deal(P(1),P(2),P(3));

503*p1 + 50*p2 + 50*p3 == 750

ans = logical
   1

100*p1 + 1309*p2 + 200*p3 == 1850

ans = logical
   1

100*p1 + 200*p2 + 407*p3 == 1600

ans = logical
   1