You are now following this question
- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.
Adding two linear inequality constraints in Optimization toolbox
5 views (last 30 days)
Show older comments
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?
Answers (3)
Alan Weiss
on 1 Jul 2020
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
9 Comments
Hussain Ja
on 1 Jul 2020
Edited: Hussain Ja
on 1 Jul 2020
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
on 1 Jul 2020
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
on 1 Jul 2020
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 ;
Hussain Ja
on 1 Jul 2020
Edited: Hussain Ja
on 1 Jul 2020
But the optimization algorithm should still be able to handle it, right? Can it be accomodated in Aineq and bineq?
Catalytic
on 1 Jul 2020
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
on 1 Jul 2020
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
on 2 Jul 2020
Edited: Matt J
on 2 Jul 2020
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
on 2 Jul 2020
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..
Walter Roberson
on 1 Jul 2020
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.
Matt J
on 2 Jul 2020
Edited: Matt J
on 2 Jul 2020
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,
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.
30 Comments
Hussain Ja
on 2 Jul 2020
So if this function is linear, how should I pass this constraint on "Power" to the main function?
Hussain Ja
on 2 Jul 2020
Edited: Hussain Ja
on 2 Jul 2020
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?
Hussain Ja
on 2 Jul 2020
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
on 2 Jul 2020
Edited: Matt J
on 2 Jul 2020
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
on 2 Jul 2020
Edited: Hussain Ja
on 2 Jul 2020
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
on 2 Jul 2020
Edited: Matt J
on 2 Jul 2020
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
on 2 Jul 2020
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
on 2 Jul 2020
You already said you were doing that in your comment here
but that it is too slow.
Matt J
on 2 Jul 2020
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.
Walter Roberson
on 2 Jul 2020
... which describes a one-level caching, similar to the formal memoize() that I talked about earlier.
Hussain Ja
on 3 Jul 2020
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
on 3 Jul 2020
Not enough information to go on.
Use
dbstop if error
dbstop if caught error
to track down where the problem is occuring.
Hussain Ja
on 3 Jul 2020
Edited: Hussain Ja
on 3 Jul 2020
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
on 3 Jul 2020
Either you objective function or your nonlinear function are emitting a cell array instead of numeric values.
Hussain Ja
on 3 Jul 2020
Edited: Hussain Ja
on 3 Jul 2020
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
on 3 Jul 2020
Edited: Walter Roberson
on 3 Jul 2020
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
on 3 Jul 2020
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);
Walter Roberson
on 4 Jul 2020
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
on 4 Jul 2020
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
on 4 Jul 2020
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
on 4 Jul 2020
Edited: Hussain Ja
on 4 Jul 2020
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
on 5 Jul 2020
Edited: Matt J
on 5 Jul 2020
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
Hussain Ja
on 5 Jul 2020
Edited: Hussain Ja
on 5 Jul 2020
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?
Hussain Ja
on 5 Jul 2020
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
on 5 Jul 2020
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
in favor of other fixed parameter passing methods described here
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!An Error Occurred
Unable to complete the action because of changes made to the page. Reload the page to see its updated state.
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)