Equality and inequality constraint

Question

Rendra Hakim hafyan 2019-9-14

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/480349-equality-and-inequality-constraint

编辑： Rendra Hakim hafyan 2019-9-15

采纳的回答： Walter Roberson

Dear All,

Currently, I am doing an integrated biorefinery process of two chemical production.

This is my objective function

p(1) = -(21655.11+61.8*x(1)-906.31*x(2)+9.34*x(2)^2)-(5924.55 -2.16*x(1)...

+ 142.75*x(3)+ 0.42*x(1)*x(3)- 0.0648*x(1)^2-0.911*x(3)^2); %% NPV

p(2) = (278.15+2.30*x(1)-6.38*x(2)+0.017*x(1)*x(2)-0.008*x(1)^2+0.075*x(2)^2)...

-(6.62+0.036*x(1)+0.21*x(3)+0.00016*x(1)*x(3)-0.00012*x(1)^2-0.0012*x(3)^2); %%TDI

p(3) =(9906.56+1660.70*x(1)-345.99*x(2)+0.941*x(1)*x(2)+0.1776*x(1)^2+3.13*x(2)^2)...

+(44362.18+25.12*x(1)^2+15.214466275*x(3)^2); %% GWP

Constraint

function [c,c_eq] = constraints(x)

c = [14.22-0.54*x(2)+0.003*x(1)*x(2)+0.005*x(2)^2-10;...

-40.52-0.022*x(1)+1.03*x(3)+0.0032*x(1)*x(3)-0.0067*x(3)^2-15]; %% This is a constraint of demand

c_eq = [(((14.22-0.54*x(2)+0.003*x(1)*x(2)+0.005*x(2)^2)/x(2))*100)-(((40.52-0.022*x(1)+1.03*x(3)+0.0032*x(1)*x(3)-0.0067*x(3)^2)/x(3))*100)- (0.3*x(1))]; %% Then this is constraint of Total glucose consumption

end

I encounter a difficulty to run my model as it only shows single dot solution. would you like to check whether there is a mistake on my model?

13 个评论
显示 11更早的评论隐藏 11更早的评论

Rendra Hakim hafyan 2019-9-14

在 MATLAB Online 中打开

Okay let me write again my problem. actually this is an integrated biorefinery of two processes

this is my new objective function

function K = biorefinery(X);
K(1) = -(-8381.33+8.599*X(1)+21635*X(2)+X(1)*X(2)*59.4+X(1)^2*0.0472+X(2)^2*-13000)...
    + (-2347.944+-4.356*X(1)+5484.462*X(3)+X(1)*X(3)*43.452+X(1)^2*-0.059+X(3)^2*-3666.667);
K(2) = (-1.578+0.947*X(1)+7.933*X(2)+X(1)*X(2)*0.16+X(1)^2*0.000133+X(2)^2*-6.667)...
    + (-502.703+1.588*X(1)+1245.998*X(3)+X(1)*X(3)*0.435+X(1)^2*0.0058+X(3)^2*-727.778);
K(3) = (40.838+2.317*X(1)+15.837*X(2)+X(1)*X(2)*0.13+X(1)^2*-0.0041+X(2)^2*-8)...
    + (-22.854+0.322*X(1)+80.064*X(3)+X(1)*X(3)*0.183+X(1)^2*-0.00077+X(3)^2*-46.667); 
end

I have got an equality and inequality constraints where for inequality constraints is demand and equality is the number of glucose.

Glucosee: (L1 + L2) = 30% x 40% x X(1)

L1 = -2.368+0.0043*X(1)+6.637*X(2)+X(1)*X(2)*0.214+X(1)^2*-0.00000267+X(2)^2*-4.667
L2 = 0.0062*X(1)+105.272*X(3)+X(1)*X(3)*0.312+X(1)^2*-0.000016+X(3)^2*-66.667

function [c,c_eq] = constraints(X)
 c = [-2.368+0.0043*X(1)+6.637*X(2)+X(1)*X(2)*0.214+X(1)^2*-0.00000267+X(2)^2*-4.667-8;...
    -41.589 -0.0062*X(1)+105.272*X(3)+X(1)*X(3)*0.312+X(1)^2*-0.000016+X(3)^2*-66.667-20];
   c_eq = (-2.368+0.0043*X(1)+6.637*X(2)+X(1)*X(2)*0.214+X(1)^2*-0.00000267+X(2)^2*-4.667...
       + -41.589 -0.0062*X(1)+105.272*X(3)+X(1)*X(3)*0.312+X(1)^2*-0.000016+X(3)^2*-66.667)...
       - 0.3*0.4*X(1);
end

how I run the model

clc
clear  
close 
%%-----------------------------------------------------------------------%%
% Optimize with gamultiobj
options = optimoptions('gamultiobj','Display','iter',...
          'MaxGeneration',1000,...
          'PopulationSize',50,...
          'CrossoverFraction',0.8,...
          'MigrationFraction',0.01,...
          'PlotFcn',@gaplotpareto);
fitness =@biorefinery;
nvars = 3;
LB = [50 0.67 0.76];
UB = [100 0.77 0.82];
ConsFcn =@constraints; 
[x,fval] = gamultiobj(fitness,nvars,[],[],[],[],LB,UB,ConsFcn,options);
% Plot results
pareto front
figure(1); 
scatter3(fval(:,1),fval(:,2),fval(:,3),'o');
xlabel('NPV ($Million)');
ylabel('FEDI');
zlabel('GWP (kg CO2-eq)');
% view(15,40)

John D'Errico 2019-9-14

编辑：John D'Errico 2019-9-14

Remember that many of the people who might be willing or able to help you are now asleep here in North America, and it is late at night for those in Europe too. Worse, this is on a weekend, so some people might not be religiously checking in. Asking for immediate help is far too demanding.

I looked quickly at your code however, and I don't understand your question. What is a single dot solution? What does that mean to you?

It looks like you have what is basically a set of three objectives, each of which seems to be quadratic in the three unknowns. You have upper and lower bounds on the variables, however the bounds seem to be surprisingly tight. So I might wonder if any solution exists at all, since you also have two inequality constraints, and an equality constraint.

You are using an optimizer, GAMULTIOBJ, so that seems logical. Your objective is composed of three sub-objectives. They are NPV, which stands for Net Present Value, perhaps? TDI, which indicates what? And what is GWP?

I might wonder if you really want to minimize net present value. But you did put a minus sign in front of that term, so maybe you handled that aspect properly.

And what do the variables X(1), X(2), and X(3) indicate?

Asking someone to check your work on a problem where we are not even given any explanation of what the variables mean? Sigh. It is you who know where those equations came from, and how they were derived, and what the variables mean.

Relax. Be patient. But it would help if you were to be more forthcoming to answer some of the questions I have brought up.

https://www.mathworks.com/matlabcentral/answers/29922-why-your-question-is-not-urgent-or-an-emergency

Rendra Hakim hafyan 2019-9-14

Dear john

I didn't mean to be pushy. just wanna get an advice only for people here where I believe they are coming from any countries not north america or the place where it's been night now only. I hope that people may be from the same region like me or country can answer it.

I am really sorry about this

okay, the result I obtained from pareto chart is only single dot where I believe that it's not the answer because for the multi-objective optimization I would have many solutions rather than one. there was something wrong on how I made my model especially for the constraints part.

actually, I already gave the new problem one

yeah K(1) is for net present value (NPV)

K(2) is global warming potential (GWP)

K(3) is fire explosion damage index (FEDI)

Since the optimization is to minimize, then my aim is to maximize the NPV so I put (-) on my NPB objective.

X(1) is variable 1 as biomass

X(2) is variable 2 as yield of product one

X(3) is variable 3 as yield of product two

I am not in rush to have the respond, I just wanna ask for those who have ever the same problem they can share the idea to me here. That's all

I don't force anyone to answer soon as well

Again, sorry for the incovenience

John D'Errico 2019-9-14

I think Walter has a good point here. That it only executes one function evaluation, because the constraint function fails due to a syntax error. Fix that, and then try again.

Rendra Hakim hafyan 2019-9-15

Dear walter and john,

Thank you for the answer

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Walter Roberson 2019-9-15

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/480349-equality-and-inequality-constraint#answer_391891

在 MATLAB Online 中打开

When you have an equality constraint, it is common to be able to get further by solving the equality for one of the variables and substituting that definition for the variable into the other portions of the function.

If you solve ceq for any one of X(1) or X(2) or X(3), you get two solutions -- that is, it is quadratic in each of the variables. We arbitrarily choose to solve ceq for the third variable X(3). We then substitute each of the two solutions in to c, getting one nonlinear inequality constraint for each of the two roots of X(3). We also substitute each of the two roots in to the fitness function, getting one fitness function for each of the two roots of X(3). When then run gamultiobj for each of the two, and then combine the results in a single graph.

Because there are upper and lower bounds on X(3), any pareto front location we find induced by the two situations might or might not meet the bounds criteria when be back-substitute the X(1) and X(2) locations into the derived equations for X(3) and check the result against the bounds. But rather than just discarding the locations that are out of bounds, we can plot them in a different color, to show the pareto front locations that were located and worked, and also the pareto front locations that were located and which did not work. We plot the locations that worked as blue downward triangles for the lower root of X(3) and as blue upward triangles for the upper root of X(3). We plot the locations that fail the X(3) bounds check as red down and upwards triangles for the two roots of X(3).

When you look at the results... you will see only red. This means that the pareto front locations for the first two variables under the equality constraint on X(3), are at places outside the bounds for X(3).

There is no solution for the combined constraints.

You can take this code and modify it to solve for a different variable in hopes that the pareto front gets lucky to find a solution, but you also need to consider the possibility that there are no pareto front locations within the constraints you posted.

I would recommend re-checking the constraint equations.

% clc
% clear  
% close 
%%-----------------------------------------------------------------------%%
% Optimize with gamultiobj
options = optimoptions('gamultiobj','Display','iter',...
          'MaxGeneration',1000,...
          'PopulationSize',50,...
          'CrossoverFraction',0.8,...
          'MigrationFraction',0.01,...
          'PlotFcn',@gaplotpareto);
fitness =@biorefinery;
nvars = 3;
LB = [50 0.67 0.76];
UB = [100 0.77 0.82];
ConsFcn =@constraints; 
[x,fval] = gamultiobj(fitness,nvars,[],[],[],[],LB,UB,ConsFcn,options);
% Plot results
pareto front
figure(1); 
scatter3(fval(:,1),fval(:,2),fval(:,3),'o');
xlabel('NPV ($Million)');
ylabel('FEDI');
zlabel('GWP (kg CO2-eq)');
% view(15,40)
X = sym('X', [1 3]);
[sc, sceq] = ConsFcn(X);
sX3 = solve(sceq, X(3));
scX3a = simplify(subs(sc, X(3), sX3(1)));   %quadratic, select smaller root
scX3b = simplify(subs(sc, X(3), sX3(2)));   %quadratic, select larger root
fscX3a = matlabFunction(scX3a, 'Vars', {X(1:2)});
fscX3b = matlabFunction(scX3b, 'Vars', {X(1:2)});
fitness3a = matlabFunction( subs(fitness(X), X(3), sX3(1)), 'Vars', {X(1:2)});
fitness3b = matlabFunction( subs(fitness(X), X(3), sX3(2)), 'Vars', {X(1:2)});
[x3a, fval3a] = gamultiobj(fitness3a, 2, [], [], [], [], LB(1:2), UB(1:2), @(X) deal(fscX3a(X),[]), options);
[x3b, fval3b] = gamultiobj(fitness3b, 2, [], [], [], [], LB(1:2), UB(1:2), @(X) deal(fscX3b(X),[]), options);
bc3a = double( subs(sX3(1), {X(1), X(2)}, {x3a(:,1), x3a(:,2)}) );
bc3b = double( subs(sX3(2), {X(1), X(2)}, {x3b(:,1), x3b(:,2)}) );
ib3a = bc3a >= LB(3) & bc3a <= UB(3);
ib3b = bc3b >= LB(3) & bc3b <= UB(3);
figure(2)
ps = 20;
ha1 = scatter3(fval3a(ib3a,1),fval3a(ib3a,2),fval3a(ib3a,3), ps, 'b', 'v');
hold on
ha2 = scatter3(fval3a(~ib3a,1),fval3a(~ib3a,2),fval3a(~ib3a,3), ps, 'r', 'v');
hb1 = scatter3(fval3b(ib3b,1),fval3b(ib3b,2),fval3b(ib3b,3), ps, 'b', '^');
hb2 = scatter3(fval3b(~ib3b,1),fval3b(~ib3b,2),fval3b(~ib3b,3), ps, 'r', '^');
hold off
xlabel('NPV ($Million)');
ylabel('FEDI');
zlabel('GWP (kg CO2-eq)');
legend([ha1, ha2, hb1, hb2], {'lower root in bounds', 'lower root out of bounds', 'upper root in bounds', 'upper root out of bounds'})