GAMULTIOBJ constraints with ODE45.

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

Answer 1

[Alan Weiss]

I believe that you are not passing the calculated values of your objective functions F correctly to your nonlinear constraint function, and that, instead, the nonlinear constraint function is keeping the control variables z positive. For an example of how to pass the calculated objective function values to the nonlinear constraint function, see Objective and Nonlinear Constraints in the Same Function.

Alan Weiss

MATLAB mathematical toolbox documentation

[Steven Taggart]

Thanks again Alan! :)

[Steven Taggart]

Alan,

I managed to sort out the problem, you were right, the optimiser was constraining the input (z) vector instead of the objectives (F) vector. I nested them within another function and only returned the necessary values to the optimiser and this has achieved what I was looking for.

Many thanks!

Steven.