Adding two linear inequality constraints in Optimization toolbox

I have the following code for maximizing revenue for a plant:

function varargout = EbsOptimize5_3(para)
para = reshape(para, 2, []);

% para = [x11,x12,x13;x21,x22,x23];
%para = reshape(para, 2, []);
global asmInfo app oc model objects Sun SunDNI SunDNImeasm M1 M2 M1measm M2measm TotalmassFlow;

if (~(size(asmInfo)))
    disp('call EbsOpen.dll')

asmInfo = NET.addAssembly('C:\Program Files\Ebsilon\EBSILONProfessional 14 P3\EbsOpen.dll');
app = EbsOpen.ApplicationClass;
oc = app.ObjectCaster;               
model = app.Open('C:\Users\Hassan Bukhari\Desktop\Priority Study\New EXP CSP 3_with_TimeSeries.ebs');
model.ActivateProfile('Charging');
objects = model.Objects;
%%%%%%%%%%%%
Sun = oc.CastToComp117(objects.Item('Sun'));
SunDNI = Sun.DNI;
%SunDNImeasm = Sun.DNI.Value;

M1 = oc.CastToComp33(objects.Item('Start_value_4'));
M2 = oc.CastToComp33(objects.Item('Start_value_5'));
M1measm = M1.M;
M2measm = M2.M;

end

SunDNImeasm = [600; 700];
powerPrice = [100; 150];
Power = zeros(1,length(SunDNImeasm));
Revenue = zeros(1,length(SunDNImeasm));
M1v = zeros(1,length(SunDNImeasm));
M2v = zeros(1,length(SunDNImeasm));
massFlow = zeros(1,length(SunDNImeasm));

for i=1:length(SunDNImeasm)
    Sun.DNI.Value =  SunDNImeasm(i);
    
    M1measm.Value = para(1,i);
    M2measm.Value = para(2,i);
    
    errors = app.NewCalculationErrors();
    model.Simulate(errors);
    
    Gen = oc.CastToComp11(objects.Item('Generator'));
    %Power = (Gen.QREAL.Value);
    Power(i) = Gen.QREAL.Value/1000 ;
    M1v(i) = M1measm.Value;
    M2v(i) = M2measm.Value;
    massFlow(i) = 3600*(M1v(i)+M2v(i));
    Revenue(i) = -(Power(i)*powerPrice(i));
end
TotalRevenue = sum(Revenue);
TotalmassFlow = sum(massFlow);
varargout{1} = TotalRevenue;
  if nargout > 1
      varargout{2} = Power;
      varargout{3} = Revenue;
      varargout{4} = TotalmassFlow;
      varargout{5} = M1v;
      varargout{6} = M2v;
  end
%TotalPower = sum(Power);
%TotalRevenue = Power*powerPrice;
%TotalRevenue = sum(Revenue);% Aineq=-ones(1,nvars);
% bineq=-9.36e3/3600;

end

I have added a linear inquality constraint on total mass flow as:

If I want to add another constraint on "Power" being calculated by the black box model attached to MATLAB, how can I do that?

Answer 1

[Alan Weiss]

You say that your new constraints are linear. In that case, you add one row to A and to b for each new constraint. For example, if there are two new constraints, x(1) + 2*x(2) <= 17 and x(3) - 4*x(4) <= -3, you would have the following:

A = -ones(1,4);
A = [A;1 2 0 0];
A = [A;0 0 1 -4];
b = -9.36e6/3600;
b = [b;17;-3];

Alan Weiss

MATLAB mathematical toolbox documentation

[Hussain Ja]

I want to place a constraint on "Power" which is being calculated by the black box model connected to MATLAB. I have two time periods, "Power" is calculated for each period. So how can I account for that? I guess I would need a constraint function. What I tried doing was copied all the code from main function and put in a new function with the following additional lines:

function [c,ceq] = csp_advisor(para)
para = reshape(para, 2, []);
c = Power(i)-170000;
ceq=[];
end

I am passing the above function as a nonlinear constraint in the Optimization tool box. It's working but its too slow to converge. Is the above code correct?

[Alan Weiss]

I do not understand the line

c = Power(i)-170000;

What is i? Where is Power calculated? What you should have is a return vector c that has two components, and then your solver tries to get both values, c(1) and c(2), to be negative.

Alan Weiss

MATLAB mathematical toolbox documentation

[Hussain Ja]

Power is calculated in the main function that I have copied above which is connected with an external black box as follows:

 model.Simulate(errors);
    
    Gen = oc.CastToComp11(objects.Item('Generator'));
    %Power = (Gen.QREAL.Value);
    Power(i) = Gen.QREAL.Value/1000 ;

[Walter Roberson]

Yup. Definitely not a linear constraint.

[Hussain Ja]

But the optimization algorithm should still be able to handle it, right? Can it be accomodated in Aineq and bineq?

[Catalytic]

How can we say anything about a black box constraint and how it may be handled? By definition, anything inside a black box cannot be examined.

[Walter Roberson]

But the optimization algorithm should still be able to handle it, right? Can it be accomodated in Aineq and bineq?

No, no chance of that.

The linear inequality and equality parameters are for dividing up the parameter space into permitted and forbidden regions through linear geometries, such as requiring that x(1) <= -2*x(2) + 5 which describes a straight line.

How can we say anything about a black box constraint and how it may be handled?

Constraints based on values that are not linear combinations of variables must be handled in the nonlinear constraint function.

As I already explained this does mean that you need to talk to the black box from both the model and the nonlinear constraint function -- and my posting explained how to do that efficiently, and explained why some of the common approaches to try to solve the difficulty are Wrnog !

[Matt J]

Constraints based on values that are not linear combinations of variables must be handled in the nonlinear constraint function.

But I think what Catalytic meant is, can we really know that to be the case? Absolutely no detail about how Power is computed is visible to us, so we do not know whether the calculations are linear combinations or not.

One question that I would add as well is, how fast is the Power computation? Are we seeing slow convergence, or is it simply that the black box calculations are slow and burdensome.

[Catalytic]

But I think what Catalytic meant is, can we really know that to be the case? Absolutely no detail about how Power is computed is visible to us, so we do not know whether the calculations are linear combinations or not.

Yes, exactly. If it turns out that Power is a linear function of the para(i), we could conceivably replace the black box calculation of the constraints with matrices..

Answer 2

[Walter Roberson]

I refer you to https://www.mathworks.com/matlabcentral/answers/557419-constraint-function-in-optimization-toolbox#answer_459217 where I explained what you need to do.

Your power is not a linear constraint: your power is being calculated based upon values created by your external process, and to put a constraint upon it you need your nonlinear constraint function to call the same external process and go through the same power calculation. I discussed the efficient way to do that in my Answer, and discussed why you need to do things that way.

Answer 3

[Matt J]

But the optimization algorithm should still be able to handle it, right? Can it be accomodated in Aineq and bineq?

One way you can test if a black box function is linear is using my func2mat File Exchange submission,

https://www.mathworks.com/matlabcentral/fileexchange/44669-func2mat-convert-linear-function-to-matrix

It will generate a matrix representation A*para(:) of the function Power(para) under the hypothesis that the function is linear. However, you would then need to verify the hypothesis by testing whether A*para(:) is in agreement with the black box calculation on lots of different choices of para.

[Hussain Ja]

So if this function is linear, how should I pass this constraint on "Power" to the main function?

[Matt J]

Using Aineq, bineq as before. Do you have reason to believe that it might be linear?

[Hussain Ja]

The relationship between "massFlow" and "Power" is inversely proportional so it's probably unlikely that the relationship is linear. However I would like to try passing it as both linear and nonlinear function.

So what specifically should I write in "Aineq"? Will "Power" have a subscript (i)?

Is there a way to pass it as a nonlinear constraint without using a new function?

[Matt J]

If the black box function PowerBlackBox(para) is equivalent to a matrix multiplication Aineq*para(:), then the constraint PowerBlackBox(para)<=bineq can be replaced by Aineq*para(:)<=bineq. That's all there is to it.

[Hussain Ja]

Shouldn't it be Aineq*Power(:)<=bineq

Sorry I am bit struggling to program this. I already have a constraint on massFlow which is:

Aineq = -ones(1,4);
bineq = -9.36e6/3600;

Would the additional constraint on "Power" be as follows:

Aineq = [Aineq;Power];
bineq = [bineq;-172000];

[Matt J]

No, the Aineq, bineq arguments are only for specifying bounds on linear combinations of the para(i). You have to be able to cast a constraint into that form in order to include it in Aineq, bineq. We know that massFlow is given by ones(1,4)*para(:), so it is indeed a linear combination of the para(i).

The question is whether there is also a different row vector v such that Power=v*para(:). If so, and if 172000 is supposed to be a lower bound on the Power, then the combined constraints would be

Aineq=[-ones(1,4) ; -v];
bineq=[-9.36e6/3600 ; -172000];

You mentioned earlier that massFlow was inversely proportional to Power, so you probably will not find such a v. However, if you actually meant that,

Power=K/massFlow

for a known constant K, then note that your lower bound on the power can be re-expressed as an upper bound on the massFlow

massFlow<=K/172000

which we already know is a linear constraint on para, and therefore you could have

Aineq=[-ones(1,4) ; ones(1,4)];
bineq=[-9.36e6/3600 ; K/172000];

[Hussain Ja]

Unfortunately, K is not known.

The problem is that there are 2 time periods in the main function with different power tarrifs i.e. $100 and $150 respectively. The first two parameters (para1 & para2) are the flows in the 1st time period. So they will generate a certain power i.e., Power(1). The last two parameters (para3 & para4) are the flows in the 2nd time period. They will generate different power, Power(2). (Power is being calculated by the black box based on the flows, so Power is an output of the black box which changes with saltFlow). I want both Power(1) and Power(2) each to be less than 172000 kW.

[Matt J]

Unfortunately, K is not known.

If the relation between massflow and power truly is

Power=K/massFlow

where K is a constant independent of para, then you could run a curvefitting routine to get K. Just run the back box to get a good range of samples of massflow and power and apply lsqcurvefit() or something.

[Hussain Ja]

I suspect the curve fitting would not be accurate as the relationship between mass flow and power changes from inverse to direct after a certain tank capacity is reached. So if we don't have a good value for K, what could an alternative approach be?

Can't we use something like

[vargout, c, ceq] = EbsOptimize5_3(para)

to pass constraint on Power?

[Matt J]

You already said you were doing that in your comment here

https://www.mathworks.com/matlabcentral/answers/557599-adding-two-linear-inequality-constraints-in-optimization-toolbox#comment_920449

but that it is too slow.

[Matt J]

Hussain Ja's comment moved here:

Yes, but I was using a separate function specifically for constraint.

Is it possible to combine the constraint function with the main function to make it computationally less expensive?

Also for some reason the constraint function that I described above works with Pattern search algorithm but not with GA.

[Matt J]

You would have to use the technique described here:

https://www.mathworks.com/help/optim/ug/objective-and-nonlinear-constraints-in-the-same-function.html

[Walter Roberson]

... which describes a one-level caching, similar to the formal memoize() that I talked about earlier.

[Hussain Ja]

I tried using the following function to combine constraint with the main function:

function [varargout,c,ceq] = EbsOptimize5_3(para)

Towards the end, I wrote the following:

c = Power - 172000;
ceq = [];

But I get the following error message:

Error running optimization.
Index in position 1 exceeds array bounds.

[Walter Roberson]

Not enough information to go on.

Use

dbstop if error
dbstop if caught error

to track down where the problem is occuring.

[Hussain Ja]

I ran the above dbstop code and got the following error message:

Caught-error breakpoint was hit in globaloptim\private\fcnvectorizer at line 13. The error was:

Conversion to double from cell is not possible.

13              y(i,:) = feval(fun,(pop(i,:)));

[Walter Roberson]

Either you objective function or your nonlinear function are emitting a cell array instead of numeric values.

[Hussain Ja]

So I changed the function to:

function [varargout,c1,c2,ceq] = EbsOptimize5_3(para)

and added the following lines:

c1 = Power(1)-170000;
c2 = Power(2)-170000;
ceq=[];

But I am still getting the same error

Caught-error breakpoint was hit in globaloptim\private\fcnvectorizer at line 13. The error was:
Conversion to double from cell is not possible.
13              y(i,:) = feval(fun,(pop(i,:)));

[Walter Roberson]

Do not use varargout as the first output parameter: it should only ever be the last parameter. varargout used like that could show up as a cell array on output.

Do not code multiple constraints by using seperate variables: all of the nonlinear inequality constraints must go into the first variable, and all of the nonlinear equality constraints must go into the second variable.

You should just not use varargout there.

EbsOptimize5_3 is your objective function. Why are you wanting to return the constraint values from it? Your objective function should return a single variable, unless you are using fmincon and the SpecifyObjectiveGradient option has been given in the options structure, in which case your objective function should return the gradient value as the second output. If you are using fmincon and you have specified HessianFcn as 'objective' and you are using trust-region-reflective then the hessian should be returned as the third output.

You are using ga or patternsearch, so you should only return a single variable from the objective function, and it must be a scalar.

Your constraint function csp_advisor should return exactly two variables, c and ceq.

If you want to return extra variables out of your functions for debugging purposes, and you are using not using fmincon with the options discussed above, then your objective function could be written something like

function [cost, power_out] = EbsOptimize5_3(para)
  ...
  if nargout > 1
    power_out = Power;
  end
end

This would not add any computation abilities, but would permit you to test your function and see what Power gets computed, by invoking EbsOptimize5_3 yourself with input para and requesting at least two outputs.

[Hussain Ja]

I tried the above suggestions and now I am getting the following error message:

Caught-error breakpoint was hit in gacreationlinearfeasible>lhsLambda at line 210. The error was:

Index in position 1 exceeds array bounds.

210     [lambda(i,:),f,e] = fmincon(fun,lambda(i,:),[],[],Aeq,beq,lb,ub,[],opts);

[Matt J]

You should make sure your objective and constraint functions work before feeding them to ga(). Call them on several test points, first. Make sure that they not only don't throw any error messages, but that they return reasonable values as well.

[Walter Roberson]

I am having difficult figuring out how that error could show up... unless possibly ga() was passed an nvars that was not a scalar integer ?

If we had the complete code we might be able to see something even without having the black box on hand, since the error looks to be occuring before the black box would be invoked.

[Hussain Ja]

My complete main code is at the top of this page.

For constraint function I am just changing the first line as follows:

function [c,ceq] = csp_advisor(para)

...............

and adding the following 2 lines at the end

c = Power(i)-170000;
ceq=[];

If I don't add all the code between function [c,ceq] = csp_advisor(para) and c = Power(i)-170000; I get errors. But when I include every line of code in between which is in the main function, then the program runs but takes too long to converge if I use Patternserarch.

For GA it gives the following error:

Caught-error breakpoint was hit in globaloptim\private\gadsplot at line 86. The error was:
Dot indexing is not supported for variables of this type.
86         position =  plotState.(solver).position; 

[Walter Roberson]

Your constraint function as posted in https://www.mathworks.com/matlabcentral/answers/557599-adding-two-linear-inequality-constraints-in-optimization-toolbox#comment_920449 (the only place you posted it) was

function [c,ceq] = csp_advisor(para)
para = reshape(para, 2, []);
c = Power(i)-170000;
ceq=[];
end

and you indicate here that you added two lines towards the bottom. That tells us that your current constraint function is

function [c,ceq] = csp_advisor(para)
para = reshape(para, 2, []);
c = Power(i)-170000;
ceq=[];
c = Power(i)-170000;
ceq=[];
end

So, now, with that constraint function csp_advisor, we know that Power is not a shared variable because your constraint function is not a nested function. So Power must be a function -- with that syntax, Power(i), that is the only two possibilities, that Power is a shared variable or a function. (I include class constructor as a kind of function.) We know that Power is not a variable being passed in to the function because we can see that the function header only passes in Para. We can see there is no global in the constraint function. So Power must be a function -- if it were not then you would get an error there because we can see it is not a variable.

Having established that Power is a function, we need the code for the function or class constructor Power.

Every time that you fail to post current code when we ask for it, we lose interest. It is frustrating for us to have to prove to you that the code you are executing cannot be the same as the problems you are describing. I am losing interest. I can tell that Matt is losing interest. Make our life easier by attaching the current source that is reproducing the problems you are talking about at this time.

[Hussain Ja]

Sorry if I didn't make myself clear. I wrote that the constraint has a few additional lines that I described above "in addition to the main function". The rest is exactly the same. So the complete constraint function is below:

function [c,ceq] = csp_advisor(para)
para = reshape(para, 2, []);

% para = [x11,x12,x13;x21,x22,x23];
%para = reshape(para, 2, []);
 global asmInfo app oc model objects Sun SunDNI SunDNImeasm M1 M2 M1measm M2measm TotalmassFlow Power;

if (~(size(asmInfo)))
    disp('call EbsOpen.dll')

asmInfo = NET.addAssembly('C:\Program Files\Ebsilon\EBSILONProfessional 14 P3\EbsOpen.dll');
app = EbsOpen.ApplicationClass;
oc = app.ObjectCaster;               
model = app.Open('C:\Users\Hassan Bukhari\Desktop\Priority Study\New EXP CSP 3_with_TimeSeries.ebs');
model.ActivateProfile('Charging');
objects = model.Objects;
%%%%%%%%%%%%
Sun = oc.CastToComp117(objects.Item('Sun'));
SunDNI = Sun.DNI;
%SunDNImeasm = Sun.DNI.Value;

M1 = oc.CastToComp33(objects.Item('Start_value_4'));
M2 = oc.CastToComp33(objects.Item('Start_value_5'));
M1measm = M1.M;
M2measm = M2.M;

end

SunDNImeasm = [600; 700];
powerPrice = [100; 150];
Power = zeros(1,length(SunDNImeasm));
Revenue = zeros(1,length(SunDNImeasm));
M1v = zeros(1,length(SunDNImeasm));
M2v = zeros(1,length(SunDNImeasm));
massFlow = zeros(1,length(SunDNImeasm));

for i=1:length(SunDNImeasm)
    Sun.DNI.Value =  SunDNImeasm(i);
    
    M1measm.Value = para(1,i);
    M2measm.Value = para(2,i);
    
    errors = app.NewCalculationErrors();
    model.Simulate(errors);
    
    Gen = oc.CastToComp11(objects.Item('Generator'));
    %Power = (Gen.QREAL.Value);
    Power(i) = Gen.QREAL.Value/1000 ;
    M1v(i) = M1measm.Value;
    M2v(i) = M2measm.Value;
    massFlow(i) = 3600*(M1v(i)+M2v(i));
    Revenue(i) = -(Power(i)*powerPrice(i));
end
TotalRevenue = sum(Revenue);
TotalmassFlow = sum(massFlow);
varargout{1} = TotalRevenue;
  if nargout > 1
      varargout{2} = Power;
      varargout{3} = Revenue;
      varargout{4} = TotalmassFlow;
      varargout{5} = M1v;
      varargout{6} = M2v;
  end
TotalPower = sum(Power);
TotalRevenue = Power*powerPrice;
TotalRevenue = sum(Revenue);% Aineq=-ones(1,nvars);

c = Power(i)-170000;
%c(2) = Power(2)-170000;
ceq=[];
end

[Matt J]

It looks like this line, which is not inside any for-loop over 'i'

c = Power(i)-170000;

should really be this

c = Power-170000;

Aside from that, though, I am no longer sure what problem this new version of your constraints is trying to address. You are evaluating almost the same code in your constraint functions as in your objective function, so you don't seem to be seeking the benefit of the caching technique I linked for you at

https://www.mathworks.com/help/optim/ug/objective-and-nonlinear-constraints-in-the-same-function.html

[Hussain Ja]

Thanks a lot. It worked!

So I tried the same constraint function by removing all the code from the main function, but that didn't speed up the results. My updated constraint function is:

function [c,ceq] = csp_advisor(para)
para = reshape(para, 2, []);
global asmInfo app oc model objects Sun SunDNI SunDNImeasm M1 M2 M1measm M2measm TotalmassFlow Power;
c = Power-170000;
ceq=[];
end

In order to get the benefit of caching technique as mentioned in the article you shared, how can I divide my objective function into two parts since my objective function is just Revenue=Power*Tarrif Rate, where Power is being calculated by the black box?

Can I write the constraint function in the same m file as the main function? Will it help speed up the results?

[Matt J]

You would put all the things that would be expensive to recompute in the computeall function. I assume that would be abll the black box calculations.

[Hussain Ja]

Why is it that I some times get the same results over and over again even when changing the input conditions? I have to quit MATLAB, re-open it and run the problem again in order to get different results, although I am doing clc; clear all in between.

[Matt J]

Who can say,without sitting next to you to see what you are running? The use of global variables is prone to causing accidents like that, though. It's highly discouraged

https://blogs.mathworks.com/videos/2010/03/08/top-10-matlab-code-practices-that-make-me-cry/

in favor of other fixed parameter passing methods described here

https://www.mathworks.com/help/optim/ug/passing-extra-parameters.html