z, z1, z2, z3 as solutions to nonlinear system of equations (many unknown variables)

Hello I have a question about a model I have been trying to solve. This is my first time using Matlab, so I may have made some beginner mistakes. I am trying to essentially solve a theoretical lagrangian (I want a symbolic solution), with 4 functions and 4 independant variables which rely on many unknown constants. I am aware that there is a Lagrangian function in Matlab, but as I am unfamiliar with the software I just wanted to solve my problem using solve().

When I try to get the solutions, my output gives solN = z, solS = z2, solW = z3 and solL = z4... I do not quite understand what this means. I am aware that this system of equations is difficult to solve, hence I have made numerous assumptions and also fixed a few constants for now (which I would prefer to keep unknown) to see if my procedure works at all. But I think I should be able to find solutions with all these simplifications...

I've also read in the forums that there used to be a bug where this happens, but I doubt that this is causing my problem...

Anyway here is my script:

syms N S W L bn f P q ip lamb theta rho bw g a p bs h b r C
%%assumptions trying to simplify the model
assume(N > 0)
assumeAlso(in(N,'integer'))
assume(S > 0)
assumeAlso(in(S,'integer'))
assume(W > 0)
assumeAlso(in(W,'integer'))
assume(bn,'real')
assume(f,'real')
assume(P,'real')
assume(q,'real')
assume(ip,'real')
assume(lamb,'real')
assume(theta,'real')
assume(rho,'real')
assume(bw,'real')
assume(g,'real')
assume(a,'real')
assume(p,'real')
assume(bs,'real')
assume(h,'real')
assume(b,'real')
assume(r,'real')
assume(C,'real')

%%setting some values for (some of the) constants
f = 90.1; %%nuclear cost per MWh LCOE in USD
g = 48; %%wind cost per MWh in USD
h = 59.1; %%Solarpv cost per MWh in USD
%%q is the amount of storage needed for 1MWh of N
%%p is the amount of storage needed for 1MWh of W
%%r is the amount of storage needed for 1MWh of S
theta = 0.5;
ip = 1/3;
rho = 1; %%particularily rho and lamb make this system complex if I leave them unknown... not sure if this matters
lamb = 1;

RE = (ip*N^lamb + (1 - ip)*(theta*W^rho + (1 - theta)*S^rho)^(lamb/rho))^(1/lamb);
PN = f*N^bn + P*q*N;
PW = g*W^bw + P*p*W^a;
PS = h*S^bs + P*r*S^b;
LAG = PN*N + PW*W + PS*S + L*(RE - C);

eqN = diff(LAG, N) == 0;
eqW = diff(LAG, W) == 0;
eqS = diff(LAG, S) == 0;
eqL = diff(LAG, L) == 0;

eqns = [eqN, eqW, eqS, eqL]
vars = [N S W L];
[solN, solS, solW, solL] = solve(eqns, vars,'Real',true)

Output:

Warning: Solutions are parameterized by the symbols: z, z1, z2, z3. To include
parameters and conditions in the solution, specify the 'ReturnConditions' value as
'true'. 
Warning: Solutions are valid under the following conditions: in(z3, 'real') & z3/3 +
48*z2^bw + z2*(48*bw*z2^(bw - 1) + P*a*p*z2^(a - 1)) + P*p*z2^a == 0 & z3/3 +
(591*z1^bs)/10 + z1*((591*bs*z1^(bs - 1))/10 + P*b*r*z1^(b - 1)) + P*r*z1^b == 0 &
z3/3 + (901*z^bn)/10 + z*(P*q + (901*bn*z^(bn - 1))/10) + P*q*z == 0 & 0 < z & 0 < z1
& 0 < z2 & z/3 - C + z1/3 + z2/3 == 0. To include parameters and conditions in the
solution, specify the 'ReturnConditions' value as 'true'. 
 
solN =
 
z
 
 
solS =
 
z1
 
 
solW =
 
z2
 
 
solL =
 
z3

Can somebody help with this?

Thanks!

Answer 1

[Walter Roberson]

Your equations have a continuous set of solutions rather than just one solution. The solution set includes all cases in which the following conditions are all simultaneously true:

48*z2^bw + z2*(48*bw*z2^(bw - 1) + P*a*p*z2^(a - 1)) + P*p*z2^a == 0
z3/3 + (591*z1^bs)/10 + z1*((591*bs*z1^(bs - 1))/10 + P*b*r*z1^(b - 1)) + P*r*z1^b == 0
z3/3 + (901*z^bn)/10 + z*(P*q + (901*bn*z^(bn - 1))/10) + P*q*z == 0
0 < z
0 < z1
0 < z2
z/3 - C + z1/3 + z2/3 == 0

Where z is N, z1 is S, z2 is W, and z3 is L

With the solution parameters matching your variables exactly, what this tells you is that you did not provide useful equations, and that the solutions basically restate the equations.

You can work a step at a time, solving one equation for one variable, and substituting that into the remaining equations, then on to solve another.

[Walter Roberson]

Ah, it looks like you can't get very far one step at a time...

[Walter Roberson]

You cannot get closed form solutions because of the variables being raised to powers. However, you can isolate equations whose roots need to be found based upon constants.

In each of the below equations, root(f(z),z) represents the set of values, z, such that f(z) = 0 -- the roots of the equation. These values need to be found numerically given the unknowns, as the exp() occurs in ways that does not lead to closed form solutions.

syms z
R1 = root(10*W^a*P*a*p - 20*exp(z)*P*q + 10*P*W^a*p + 480*W^bw*bw - 901*exp(z*bn)*bn - 901*exp(z*bn) + 480*W^bw, z);
R3 = root(10*W^a*P*a*p - 10*P*exp(z*b/bs)*b*r - 10*P*exp(z*b/bs)*r + 10*P*W^a*p + 480*W^bw*bw - 591*exp(z)*bs - 591*exp(z) + 480*W^bw,z);
R2 = root(10*(-exp(R1) + 3*C - exp(R3/bs))^a*P*a*p - 20*exp(z)*P*q + 10*P*(-exp(R1) + 3*C - exp(R3/bs))^a*p + 480*(-exp(R1) + 3*C - exp(R3/bs))^bw*bw - 901*exp(z*bn)*bn - 901*exp(z*bn) + 480*(-exp(R1) + 3*C - exp(R3/bs))^bw,z);
R4 = 10*(-exp(R1) + 3*C - exp(R3/bs))^a*P*a*p - 10*P*exp(z*b/bs)*b*r - 10*P*exp(z*b/bs)*r + 10*P*(-exp(R1) + 3*C - exp(R3/bs))^a*p + 480*(-exp(R1) + 3*C - exp(R3/bs))^bw*bw - 591*exp(z)*bs - 591*exp(z) + 480*(-exp(R1) + 3*C - exp(R3/bs))^bw,z);

Then

N = exp(R2);
S = exp(R4/bs);
L = -3*P*p*(a + 1)*(-exp(R1) + 3*C - exp(R3/bs))^a - 144*(-exp(R1) + 3*C - exp(R3/bs))^bw*(bw + 1);
W = -exp(R1) + 3*C - exp(R3/bs);

[Valentin Gerschbacher]

Thank you for your detailed help and fast response, unfortunately this is not the kind of answer I am looking for, so instead I will need to change the model.