Equality and inequality constraint

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?

[Walter Roberson]

How are you running your model? And what do you mean by a single dot solution?

[Rendra Hakim hafyan]

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)

[Rendra Hakim hafyan]

Is there anyone who can help me out?

I just wanna have an advice whether I have been right or wrong in solving such problem

Thank you for your kindness

[Walter Roberson]

Sorry, my partner would be pretty upset if I got out of bed at 02:30 to paste in all the code to experiment with.

[Rendra Hakim hafyan]

Sorry, I don't get it

I just wanna know which part that I make the mistake

what I do is already right or do I need to fix some parts?

may anyone else here can help me?

[Walter Roberson]

It is the middle of the night in North America and to test your code I would need to go over to my office (it is too complicated for command line remote access from my cell phone)

My working assumption is that your constraints are either wrong or very restrictive.

[Rendra Hakim hafyan]

Ah I see

sorry, I don't even realize about the time

okay never mind

I think so

[John D'Errico]

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]

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

[Walter Roberson]

In your constraints,

 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];

The first of those lines has one expression in the row. The second of those two lines has two expressions in the row. That would generate an error, preventing more than one iteration.

-41.589 -0.0062*X(1)etc

is not subtraction: it is -41.589 as one element, and -0.0062*X(1)etc as a second element. Spaces matter in [] lists. Either remove the space before -0.0062*X(1)etc or else put a space after the -, making -41.589 - 0.0062*X(1)etc

[Walter Roberson]

I turn out to be the latest-working North American volunteer, but especially after 06:00 UTC (01:00 Central time) I am increasingly likely to work from bed on my cell phone. On weekdays, the regular European and Indian volunteers start coming online about 08:00 UTC but historically they are a bit irregular until about 10:00 UTC. Between about 06:00 UTC and 09:00 UTC on weekdays the regular volunteers tend not to be around, or tend to be me working extra late (and getting too tired for complex answers) or working on my cell phone from bed.

[John D'Errico]

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]

Dear walter and john,

Thank you for the answer

Answer 1

[Walter Roberson]

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'})

[Rendra Hakim hafyan]

yeah I totally agree with you about the equality costraints.

to be honest, I would like to integrate two processes where you can see the picture below.

firstly, I make an empirical model of each objective function with function of X1,X2, and X3. also for the model of product demand (p1 and p2)

X1 : biomass
X2: yield of HMF dehydration
X3: yield of fermentation

Inequality constraint

P1 = funciton (X1,X2) <= the demand of product P1
P2 = function(X1,X3) <= the demand of product P2

Equality constraint

Since,

L1 = 0.4 X1
L2 = 0.3 L1
L2 = L3 + L4

since, L3 and L4 relate to P1 & P2 so,

P1 + P2 = 0.3 x 0.4 x X1 

Just to clear my problem up. Do you think my model has been right?

okay walter, I'd like to check again my equality constraint

[Walter Roberson]

P1 = funciton (X1,X2) <= the demand of product P1

Typically one would express that the demand on a product must be less than or equal to the rate at which it is produced ?

[Rendra Hakim hafyan]

yeah, I would like to make sure that the P1 should not exceed the demand to avoid product excess