Hello,
I have had assistance on here recently with this problem but I am still facing issues so I'd like to include some more information this time, I am fairly new to Matlab and optimisation so I am encountering some problems. I have a function that models valve dynamics using the ODE45 solver. This function takes an input vector [z] that contains x & y coordinates, these create a curve that determines how the valve behaves. The output of a function is a 2 element vector (F) that I am trying to minimise using the GAMULTIOBJ optimisation tool, the values in [F] need to be kept positive. I had some help on here to create a non-linear constraint function to keep the values positive but it does not seem to help, my Pareto front from the GAMULTIOBJ is always outside where I am trying to constrain it. I can choose values for the input vector that can produce the positive values I am looking for so the function works fine but it is the optimisation that seems to be the problem. I have included the function below. The function used by ODE45 is included for completeness.
function [F] = runoptimise(opts) FitnessFunction = @safety_relief_dynamics; ConstraintFunction = @mycon; numberofVariables = 12; lb = []; ub = []; [F] = gamultiobj(FitnessFunction, numberofVariables,[],[],[],[],lb,ub,ConstraintFunction,opts); function F = safety_relief_dynamics(z) tspan=[0,60]; %%time span
ic=[3.25e5,296,0,0,0,0,60]; %%initial conditions
options=odeset('AbsTol',1e-8,'RelTol',1e-11,'InitialStep',0.00000001); %%error limits to control time step
[t,y]=ode45(@(t,y) broady(t,y(1),y(2),y(3),y(4),y(5),y(6),y(7),z,tspan(2)),tspan,ic,options); P = y(:,1); P_max=max(P); OP = ((((P_max) - 3.3e5) / 3.3e5)*100)-1.5; BD = (((3.3e5 - (P(end)))/ 3.3e5)*100)-1.5; F = [OP, BD]; %assignin('base', 'F', F);
fprintf ('\n'); fprintf ('Blowdown : %.3f',BD); fprintf ('\n'); fprintf ('Overpressure : %.3f', OP); fprintf ('\n'); end function [c,ceq] = mycon (F) c(1)= -F(1); c(2)= -F(2); ceq = []; end end%%function used by ODE45
function y = broady(t,P,T,~,~,x,ve,~,z,tfinal) %Timer
fprintf('TIME: %.3f of %.3f\n',t,tfinal) %constants
V = 1.5; %%tank volume
g = 9.81; %%gravitational constant
R=287; %%ideal gas constant
Cp=1005; %%specific heat capacity of air
gamma=1.4; %%heat capacity ratio
d=0.01393; %%Inlet diameter
A_i=pi*(d/2)^2; %%Inlet area
P_atm=101000; %%atmospheric pressure
P_set=3.3e5; %%set pressure
T_atm=293; %%atmospheric temperature
T_o=288; %%Inlet gas temperature
ma=0.651; %%mass of moving parts
k=16000; %%spring stiffness
c=1; %%Damping coefficient
%%Discharge coefficient and flow area
CD=0.2; A=3.14*d*x; if A>A_i A=A_i; end if A <0 A = 0; end %choked/subsonic conditions
if Patm/P>0.5283 P_prime=P_atm; T_prime=T_atm; else P_prime=0.5283*P; T_prime=0.833*T; end %Mass flow rate(algebraic eqn) and conditions
m_e=A*CD*(P_prime/(R*T_prime))*(2*Cp*T*(1-(P_prime/P)^((gamma-1)/gamma)))^0.5; %%mass outflow rate
if (1-P_prime/P)<=0 m_e=0; end if t <30 m_i = 0.02; else m_i = 0; end m_v=m_i-m_e; %%mass rate of change of vessel
%odes
dP=(1/V)*(Cp*(gamma-1)*(m_i*T_o-m_e*T)); %%change in pressure
dT=(T/P)*dP-((R*T^2)/(P*V))*(m_v); %%change in temperature
dm_e=m_e; %%amount of mass to leave vessel
dm_v=m_v; %%mass change of vessel
dx=ve; %%rate of change of position
%Valve opening and closing limits
if x>0.004 if dx>0 dx=0; end end if x<1e-10 if dx<0 dx=0; end end %Force constants
x_1 = z(1:6); y_1 = z(7:12); f_1 = pchip(x_1,y_1,x); F=(f_1*P)-(y_1(1)*P_set); nf = F; dve=(F-c*ve-k*x)/(ma); %%rate of change of disc velocity
if x < 1e-10 if dve < 0 dve = 0; end end %matrix/vector
y=[dP;dT;dm_e;dm_v;dx;dve;nf]; end%%script to run GAMULTIOBJ
opts = optimoptions('gamultiobj', 'UseParallel', true, 'PlotFcn',@gaplotpareto); rng default; runoptimise(opts)
The constraint function @mycon does not seem to hold the [F] values positive, whither I have this constraint function active or not, it seems to make no difference and I get the same Pareto front from the optimiser, does the constraint look correct? Is it applied correctly? I have included a pic of the Pareto front showing it is not in the x-y positive region.
Thanks for your time! 🙂
Steven
Best Answer