Sequential Constraint in Optimization algorithm.
1 Ansicht (letzte 30 Tage)
Ältere Kommentare anzeigen
Pavitra Jain
am 24 Aug. 2022
Beantwortet: Walter Roberson
am 26 Aug. 2022
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);
2 Kommentare
Alan Weiss
am 25 Aug. 2022
What is your question?
Alan Weiss
MATLAB mathematical toolbox documentation
Akzeptierte Antwort
Walter Roberson
am 26 Aug. 2022
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".
0 Kommentare
Weitere Antworten (0)
Siehe auch
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!