Sequential Constraint in Optimization algorithm.

I am working on a problem where i have 2 constraints, 2nd contraint is a subset of 1st constraint and to calculate the output for 2nd constraint we need to satisfy the 1st constraint and then the 2nd constraint otherwise we will get error value.

For now i went for a round about way to do this by putting a if condition that bypass that calculation if the first constraint is not satisfied.

I have attached the code below:

accelero4poly is my objective function

fanconst4poly is my non linear constraints - here c value is for 1st constraint that i mentioned and ceq is 2nd constraint.

function y = accelero4poly(x)

syms X 
rho = 2330;  %Kg/m3 Density
E = 169e9; %Pa Elastic Modulus
gap = 1E-6; %metres 1e-6 Gap between Electrodes
e = 8.8581e-12; %SI Dielectric constant
th = 20e-6; %metres Thickness of beam

g = 9.8; %1g Acceleration acting on the Beam

L = x(5); %Length of the beam

w = @(X) x(1)*X.^3 + x(2)*X.^2 + x(3)*X + x(4) + x(6)*X.^4; %metres  Polynomial for Width

dwdX = @(X) 3*x(1)*X.^2 + 2*x(2)*X + x(3)+ 4*x(6)*X.^3; %unitless Polynomial for the derivative of width

%--If it exceeds then function returns a 0-----%

I0 = (1/12)*(th)*(w(0))^3;

wx = @(X) (x(1)*X^3 + x(2)*X^2 + x(3)*X + x(4)+ x(6)*X.^4)*X; %Declaring width times x

%Moment acting at the fixed end - as (d2y/dx2)(x = 0) = M/EI0 =
%(density*thickness*Acc/EI0)*(integration of width times x wrt x for 0 to L)
M = double((rho*th*g/(E*I0))*(int(wx,X,0,L)));


%Now comes shear force at the fixed end - as (d3ydx3)(x=0) = (1/EI0)*(-(density*thickness*Acc*integration of w wrt x from 0 to L) - E*(thickness/4)*(width)^2*(der of w wrt x)*d2ydx2)

V = double((1/(E*I0))*(-(rho*th*g*int(w,X,0,L))-E*(th/4)*(w(0))^2*(dwdX(0))*M));

xRange = [0,L];
yinit = [0 0 M V];
sol = ode45(@(t,Y) deflection(t,Y,x) ,xRange,yinit);
xinit = linspace(0,L,300);
Sxinit = deval(sol,xinit);
plot(xinit,Sxinit(1,:));

% Length section along the width is ds = (1+(dydx)^2)^0.5
ds = (1 + (dwdX(xinit)).^2).^0.5;

C2 = ((1E20*e*th).*ds)./(gap - Sxinit(1,:));

Caftsimps = simps(xinit,C2);

C1 = ((1E20*e*th).*ds)./(gap);

Cbefsimps = simps(xinit,C1);

errorsimps = abs(Caftsimps - Cbefsimps);

% ----Check if resonant freq lies in the required limit----%
wydx = simps(xinit,w(xinit).*Sxinit(1,:));
wy2dx = simps(xinit,w(xinit).*(Sxinit(1,:).^2));
freq = (g*wydx/wy2dx)^0.5;
display(freq)

%-----If it exceeds then function returns a 0-----%

y = -1*errorsimps;

end

function dYdt = deflection(t,Y,x)
Y0 = Y(1);
Y1 = Y(2);
Y2 = Y(3);
Y3 = Y(4);

dY0dt = Y1;
dY1dt = Y2;
dY2dt = Y3;

rho = 2330;  %Kg/m3 Density
E = 169e9; %Pa Elastic Modulus
g = 9.8; %1g Acceleration acting on the Beam

%d4y/dx4 = (1/w^2)*(-6*w*dwdx*d3ydx3 - (6*w*dwdx^2+3*w^2*d2wdx2)*d2ydx2 + 12*density*acc/E)
recw2 = (1/((x(1)*t.^3 + x(2)*t.^2 + x(3)*t + x(4)+ x(6)*t.^4)^2));
w = (x(1)*t.^3 + x(2)*t.^2 + x(3)*t + x(4)+ x(6)*t.^4);
dwdx = (3*x(1)*t.^2+2*x(2)*t+x(3)+ 4*x(6)*t.^3);
d2wdx2 = (6*x(1)*t + 2*x(2)+ 12*x(6)*t.^2);

dY3dt = recw2*(-6*w*dwdx*Y3 - (6*w*(dwdx)^2 + 3*(d2wdx2)*(w)^2)*Y2 + 12*rho*g/E);


dYdt = [dY0dt ; dY1dt ; dY2dt ; dY3dt];
end

function [c, ceq] = fabconst4poly(x)

syms X 
rho = 2330;  %Kg/m3 Density
E = 169e9; %Pa Elastic Modulus
th = 20e-6; %metres Thickness of beam

g = 9.8; %1g Acceleration acting on the Beam

L = x(5); %Length of the beam



w = @(X) x(1)*X.^3 + x(2)*X.^2 + x(3)*X + x(4) + x(6)*X.^4; %metres  Polynomial for Width

wmax = max(w(linspace(0,L)));
wmin = min(w(linspace(0,L)));
c = [(wmax - 10E-6) , (2E-6 - wmin)];

idx = find(c>0);
if idx > 0
   ceq = 2;
    return
end

dwdX = @(X) 3*x(1)*X.^2 + 2*x(2)*X + x(3)+ 4*x(6)*X.^3; %unitless Polynomial for the derivative of width

%-----If it exceeds then function returns a 0-----%

I0 = (1/12)*(th)*(w(0))^3;

wx = @(X) (x(1)*X^3 + x(2)*X^2 + x(3)*X + x(4)+ x(6)*X.^4)*X; %Declaring width times x

%Moment acting at the fixed end - as (d2y/dx2)(x = 0) = M/EI0 =
%(density*thickness*Acc/EI0)*(integration of width times x wrt x for 0 to L)
M = double((rho*th*g/(E*I0))*(int(wx,X,0,L)));


%Now comes shear force at the fixed end - as (d3ydx3)(x=0) = (1/EI0)*(-(density*thickness*Acc*integration of w wrt x from 0 to L) - E*(thickness/4)*(width)^2*(der of w wrt x)*d2ydx2)

V = double((1/(E*I0))*(-(rho*th*g*int(w,X,0,L))-E*(th/4)*(w(0))^2*(dwdX(0))*M));

xRange = [0,L];
yinit = [0 0 M V];
sol = ode45(@(t,Y) deflection(t,Y,x) ,xRange,yinit);
xinit = linspace(0,L);
Sxinit = deval(sol,xinit);
plot(xinit,Sxinit(1,:));


% ----Check if resonant freq lies in the required limit----%
wydx = simps(xinit,w(xinit).*Sxinit(1,:));
wy2dx = simps(xinit,w(xinit).*(Sxinit(1,:).^2));
freq = (g*wydx/wy2dx)^0.5;

ceq = freq - 2E6;

%-----If it exceeds then function returns a 0-----%

end

function dYdt = deflection(t,Y,x)
Y0 = Y(1);
Y1 = Y(2);
Y2 = Y(3);
Y3 = Y(4);

dY0dt = Y1;
dY1dt = Y2;
dY2dt = Y3;

rho = 2330;  %Kg/m3 Density
E = 169e9; %Pa Elastic Modulus
g = 9.8; %1g Acceleration acting on the Beam

%d4y/dx4 = (1/w^2)*(-6*w*dwdx*d3ydx3 - (6*w*dwdx^2+3*w^2*d2wdx2)*d2ydx2 + 12*density*acc/E)
recw2 = (1/((x(1).*t.^5 + x(2).*t.^4 + x(3).*t.^3 + x(4).*t.^2 + x(5).*t + x(6))^2));
w = (x(1).*t.^5 + x(2).*t.^4 + x(3).*t.^3 + x(4).*t.^2 + x(5).*t + x(6));
dwdx = (5*x(1).*t.^4 + 4*x(2).*t.^3 + 3*x(3).*t.^2 + 2*x(4).*t + x(5));
d2wdx2 = (20*x(1).*t.^3 + 12*x(2).*t.^2 + 6*x(3).*t + 2*x(4));

dY3dt = recw2*(-6*w*dwdx*Y3 - (6*w*(dwdx)^2 + 3*(d2wdx2)*(w)^2)*Y2 + 12*rho*g/E);


dYdt = [dY0dt ; dY1dt ; dY2dt ; dY3dt];
end

function [x,fval,exitflag,output,population,score] = untitled(nvars,lb,ub)
%% This is an auto generated MATLAB file from Optimization Tool.

%% Start with the default options
options = optimoptions('ga');
%% Modify options setting
options = optimoptions("ga",'PlotFcn',{@gaplotbestf,@gaplotmaxconstr}, ...
    'Display','iter',opts,'UseParallel',true);
[x,fval,exitflag,output,population,score] = ...
ga(@accelero4poly,nvars,[],[],[],[],lb,ub,@fabconst4poly,options);

ub = [900*(1E-6*1E24) , 900*(1E-6*1E18) , 900*(1E-6*1E12) , 900 , 150E-6 , 10E-6];
lb = [-900*(1E-6*1E24) , -900*(1E-6*1E18) , -900*(1E-6*1E+12) , -900 , 50E-6 , 2E-6 ];
nvars = 6;
[x , fval] = untitled(nvars,lb,ub);
display(x);
display(fval);

[Alan Weiss]

What is your question?

Alan Weiss

MATLAB mathematical toolbox documentation

[Pavitra Jain]

I want to add two constraint where it moves to 2nd constraint only when 1st constraint is satisfied.

Answer 1

[Walter Roberson]

For now i went for a round about way to do this by putting a if condition that bypass that calculation if the first constraint is not satisfied.

That is fine. As long as each time, it emits the same number of entries for c as previous times, and the same number of entries for ceq as previous times (and it is advised that any "used" constraint occupies a consistent position in the output order.)

ga() and the other optimizers do not have any built-in way to express "if all of this set of constraints were satisfied, then call this for additional constraints".