Neeta Dsouza answered . 2024-05-17 19:47:09

function [x,f,eflag,outpt] = myModelParamsSolver

xLast = []; % Last place computeall was called
myf = []; % Use for objective at xLast
myc = []; % Use for nonlinear inequality constraint
myceq = []; % Use for nonlinear equality constraint
N = 10;
x0 = rand(1,2+2*N);
fun = @objfun; % the objective function, nested below
cfun = @constr; % the constraint function, nested below
options = optimoptions('fmincon');
lb = zeros(1,2+2*N);
ub = ones(1,2+2*N);
% Call fmincon
[x,f,eflag,outpt] = fmincon(fun,x0,[],[],[],[],lb,ub,cfun,options);

      function y = objfun(x)
          if ~isequal(x,xLast) % Check if computation is necessary
              [myf,myc,myceq] = myDiscreteTimeModel(x);
              xLast = x;
          end
          % Now compute objective function
          y = myf;
      end

      function [c,ceq] = constr(x)
          if ~isequal(x,xLast) % Check if computation is necessary
              [myf,myc,myceq] = myDiscreteTimeModel(x);
              xLast = x;
          end
          % Now compute constraint functions
          c = myc; % In this case, the computation is trivial
          ceq = myceq;
      end

      function [f,c,ceq] = myDiscreteTimeModel(x)
          p = zeros(1,N);
          v = p;
          u = p;
          Ts = 0.01;
          p_max = 100;
          p_init = 10;
          vel_init = 1;

          % initial conditions:
          p(1)=p_init;
          v(1)=vel_init;
          p1 = x(1);
          p2 = x(2);
          v1 = x(3:3+N-1);
          v2 = x(3+N:end);

          for n = 1:N
              t = N*Ts;

              % non linear equalities:
              s_star= v1(n) * (v1(n)-v2(n));
              % where v1 and v2 are parameters of entities 1 and 2 respectively at nth instant.
              % * this is the non-linear equality/constraint.
              %s_star2 = v1(n) * (v1(n)-v1(n-1));
              % Note: we use either s_star or s_star2 in our code

              % linear equalities:
              p(n+1)= p(n)+ v(n)*t + 0.5*u(n)*t^2;
              v(n+1)= v(n)+ u(n)*t;
              u(n+1)= (s_star / (p1-p2))^2 ; % p1 and p2 correspond to entity 1 and 2's p parameter
          end

          % end conditions:
          ceq = v(N);
          c = p(N)-p_max;
          f = sum(u(:));
      end
  end

Not satisfied with the answer ?? ASK NOW