How to set upper and lower bounds on vpasolve

syms P1 P4 P5 lambda
Load_hour_2 = 470;
%equations
F1 = 0.0013*P1^2+11.154*P1+60;
F4 = 0.0016*P4^2+11.336*P4+52;
F5 = 0.0016*P5^2+14.541*P5+52;
%calculate partial derivatives
F1_D = diff(F1,P1);
F4_D = diff(F4,P4);
F5_D = diff(F5,P5);
%set upper and lower bounds?
assume(P1>=90 & P1<=270);
assume(P4>=50 & P4<=180);
assume(P5>=50 & P5<=180);
%the equation below gives me a soution but does not take in to consideration the assumptions made above and the solution is not viable for what i am trying to accomplish
[lambda_calc,gen1,gen4,gen5] = vpasolve(F1_D==lambda,F4_D==lambda,F5_D==lambda,(P1+P4+P5)==Load_hour_2,[P1, P4, P5,lambda]);
%setting the bounds this way results in a sym 0-by-4 result
[lambda_calc,gen1,gen4,gen5] = vpasolve(F1_D==lambda,F4_D==lambda,F5_D==lambda,(P1+P4+P5)==Load_hour_2,[P1, P4, P5,lambda],[90 270;50 180;50 180;0 100]);

I am trying to calculate the economic dispatch for my unit commitment and I cannot seem to get an answer that follows my upper and lower bounds for the solution, they are either ignored or I recevie a 'sym 0-by-4' result. Is there a different type of solve function that would work for this? (lambda does not need to be constrained it just gives errors if I dont specify anything for that variable)

Answer 1

[Walter Roberson]

vpasolve() has a ... challenge... with respect to bounds.

vpasolve() does a modified Newton-Raphson. Those kinds of methods inherently involve projection and refinement, with projection to both sides of the actual root being inherent to the process.

However... in vpasolve(), if the projection extends past the limits that were set, then vpasolve() gives up. vpasolve() does not move to the closest boundary and try again: it just gives up.

The location of projection can get pretty wild in Newton-Raphson methods. If you happen to guess near to the actual root, then the direction to the actual root can be quite steep, resulting in a large projection... that can then result in vpasolve giving up even though it was quite close to a root in a recent step.

[Walter Roberson]

By the way, one consequence of this behaviour of Newton-Raphson type calculations is that sometimes you can get faster convergence by starting further away from the correct root!

However, when you use limits for vpasolve() you cannot specify starting points. You can control whether vpasolve() initialized randomly or not, but you cannot control the locations it starts with if you specify limits.

[Mackenzie Ethier]

Okay thank you for explaining that! Would you have any recommendations for what sort of solve I could use with constraints?

[Walter Roberson]

Sure. For that problem, the solver you should use

syms P1 P4 P5 lambda
Load_hour_2 = 470;
%equations
Q = @(v) sym(v);
F1 = Q(0.0013)*P1^2 + Q(11.154)*P1 + 60;
F4 = Q(0.0016)*P4^2 + Q(11.336)*P4 + 52;
F5 = Q(0.0016)*P5^2 + Q(14.541)*P5 + 52;
%calculate partial derivatives
F1_D = diff(F1,P1);
F4_D = diff(F4,P4);
F5_D = diff(F5,P5);
%set upper and lower bounds?
%assume(P1>=90 & P1<=270);
%assume(P4>=50 & P4<=180);
%assume(P5>=50 & P5<=180);
%the equation below gives me a soution but does not take in to consideration the assumptions made above and the solution is not viable for what i am trying to accomplish
sol = solve(F1_D==lambda,F4_D==lambda,F5_D==lambda,(P1+P4+P5)==Load_hour_2);

lambda_calc = sol.lambda;
gen1 = sol.P1;
gen4 = sol.P4;
gen5 = sol.P5;

You will notice that the results violate the desired constraints. That is because the equations are linear equations and there is only one solution and that one solution is outside the constraints.

Derivative of an equation of degree 2 is an expression of degree 1. Sum of expressions of degree 1 is degree 1. So all parts are linear, and the solution is whatever it is.

Now, if this is an optimization problem, and the equations give the optimum, then the fact that the solution to the equation is not inside the constraints implies that the minimum lies at the boundary of at least one of the constraints.

My tests show that the constrained solution is at P4 = 180, P5 = 50, and is P1 = 18000279825/75000338, lambda = 959779497597/75000338000

[Mackenzie Ethier]

Thank you very much!