Sequential Constraint in Optimization algorithm.

4 次查看（过去 30 天）
Pavitra Jain
Pavitra Jain 2022-8-24
回答： Walter Roberson 2022-8-26
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 个评论
Alan Weiss
Alan Weiss 2022-8-25
What is your question?
Alan Weiss
MATLAB mathematical toolbox documentation
Pavitra Jain
Pavitra Jain 2022-8-26
I want to add two constraint where it moves to 2nd constraint only when 1st constraint is satisfied.

请先登录，再进行评论。

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

采纳的回答

Walter Roberson
Walter Roberson 2022-8-26
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 个）

类别

在 Help Center 和 File Exchange 中查找有关 Assembly 的更多信息

标签

产品


版本

R2022a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by