Hi all; I try to use fmincon in order to find the optimal control force F for an active suspension system, every time I run fmincon, F and the objective function did not change. Kindly, can any on explain to me why I get this problem. Thanks

Where the system is a dynamic system and it is subjected to the equation of motion. The objective function is determined based on Lagrange integral parameterfun(F).

Fo=0;
lb=0;
ub=3000;  
options = optimset('Display','iter','TolFun',1e-8,'algorithm','sqp');

[Fopt,fval,exitflag,output] = fmincon(@(F) parameterfun(F),Fo,...
                                 [],[],[],[],lb,ub,[],options);

where parameterfun(F) is the objective function which calculate the value of F

    function [val,F] = parameterfun(F)
    val = 0;
    [xu,t,z0dot]=xu_det(F);
    xu_frame=xu;
    input=z0dot;
    r1 = 1e+5;
    r2 = 0.5;
    q = 1e-5;
    R = diag([r1, 0, 0, 0, q, r2]);

    xu_frame_zs = zeros(6,length(t));

    for i=0:0.005:(t-0.005)
        dx = f(t(i),xu_frame(1:4,i), xu_frame(5,i),input(i));
        xu_frame_zs(:,i) = [xu_frame(:,i); dx(4)];
        val = val+ (t(i+0.005)-t(i))*xu_frame_zs(:,i)'*R*xu_frame_zs(:,i);
    end
    end

where xu_det(F) is a function which is:

      function [xu,t,z0dot] = xu_det(F)
      param.ms = 325;               % 1/4 sprung mass (kg)
      param.mus = 65;               % 1/4 unsprung mass (kg)
      param.kus = 232.5e3;          % tire stiffness (N/m)

      ks=2.36e4;
      cs=1.03e3;

      Aqcar = [0 1 0 0;-param.kus/param.mus -cs/param.mus ks/param.mus cs/param.mus;0 -1 0 1;0 cs/param.ms -ks/param.ms -cs/param.ms];
      Bqcar = [0 -1;-1/param.mus 0; 0 0; 1/param.ms 0];
      Cqcar = [1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1];
      Dqcar = 0;

      qcar = ss(Aqcar,Bqcar,Cqcar,Dqcar);

      x0 = [0 0 0 0]';                        % initial state
      load IRI_737b                           % road profile data
      ddx = road_x(2) - road_x(1);             % spacial step for input data
      v = 10;                 % vehicle velocity (m/s)
      dt = 0.005;
      dt2 = ddx/v;                             % time step for input data
      z0dot = [0 diff(road_z)/dt];            % road profile velocity     
      tmax = 5;                               % simulation time length
      t = 0:dt:tmax;
      x = v*t;                 % time/space steps to record output
      z0dot_f=zeros(length(z0dot),2);

      z0dot_f(:,1)=z0dot;
      z0dot_f(:,2)=F;
      u = interp1(road_x,z0dot_f,x);%umf = 3;
      umf = 0.01;
   % Simulate quarter car model
      y = lsim(qcar,u*umf,t,x0);   
      yu=y;
      yu(:,5)=F;
      xu=yu';
    end

and f() is another function is used to determine the states variables derivative dx at each F value

function dx = f(~, x,F,z0dot )
J = jacobian();
A = J(1:4,1:4);
B = J(1:4, 5);
b_ramp = [-1;0; 0; 0];
    dx= A* x + B * F + b_ramp * z0dot; 
end

and jacobian is another function is used to determine A and B matrices

function J = jacobian()
param.ms = 325;               % 1/4 sprung mass (kg)
param.mus = 65;               % 1/4 unsprung mass (kg)
param.kus = 232.5e3;          % tire stiffness (N/m)

ks=2.36e4;
cs=1.03e3;

A = [ 0 1 0 0;
    -param.kus/param.mus -cs/param.mus ks/param.mus cs/param.mus;
    0 -1 0 1;
    0 cs/param.ms -ks/param.ms -cs/param.ms];
B = [0 -1/param.mus 0 1/param.ms]'; 
J = [A, B];
end

when I run fmincon the results as:

                                                           Norm of First-order
   Iter F-count            f(x) Feasibility  Steplength        step  optimality
      0       2   0.000000e+000  0.000e+000                          1.000e+000
      1       4   0.000000e+000  0.000e+000  1.000e+000  0.000e+000  0.000e+000

Local minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the selected value of the function tolerance,
and constraints are satisfied to within the default value of the constraint tolerance.

<stopping criteria details>

Answer 1

[Alan Weiss]

Did you try plotting your objective function to make sure that it is not constant near the initial point? Are you sure that you calculated the gradient of your function correctly? If so, then tell fmincon that you are supplying the gradient by setting

options.GradObj = 'on';

If you are not supplying a gradient, then you might need larger finite difference steps sizes, as described in the documentation on optimizing simulations.

Good luck,

Alan Weiss

MATLAB mathematical toolbox documentation

[Muna Shehan]

Thanks for your reply. I did not plot the objective function. Kindly; can you guide me toward a matlab example for more explanation. Thanks a lot

[Alan Weiss]

It seems that you have a problem with just a one dimensional independent variable that you call F. So why not try to evaluate parameterfun(F) for various values of F and see whether or not it changes, especially for small F near your initial point 0, such as 1e-8.

Also, you might as well use fminbnd for a 1-D problem, as fmincon is sort of overkill. But suit yourself.

Alan Weiss

MATLAB mathematical toolbox documentation

[Muna Shehan]

Thanks for your reply, when i use fminbnd the results which i get are as showen:

 Func-count     x          f(x)         Procedure
    1         1145.9            0        initial
    2         1854.1            0        golden
    3         2291.8            0        golden
    4        2562.31            0        golden
    5        2729.49            0        golden
    6        2832.82            0        golden
    7        2896.67            0        golden
    8        2936.14            0        golden
    9        2960.53            0        golden
   10        2975.61            0        golden
   11        2984.93            0        golden
   12        2990.68            0        golden
   13        2994.24            0        golden
   14        2996.44            0        golden
   15         2997.8            0        golden
   16        2998.64            0        golden
   17        2999.16            0        golden
   18        2999.48            0        golden
   19        2999.68            0        golden
   20         2999.8            0        golden
   21        2999.88            0        golden
   22        2999.92            0        golden
   23        2999.95            0        golden
   24        2999.97            0        golden
   25        2999.98            0        golden
   26        2999.99            0        golden
   27        2999.99            0        golden
   28           3000            0        golden
   29           3000            0        golden
   30           3000            0        golden
   31           3000            0        golden
   32           3000            0        golden
   33           3000            0        golden
   34           3000            0        golden
   35           3000            0        golden
   36           3000            0        golden
   37           3000            0        golden

Optimization terminated: the current x satisfies the termination criteria using OPTIONS.TolX of 1.000000e-012

>> in other words i still have a problem, because even if i change ub to a different value the optimizer run until the ub value and return ub as optimal value

[Alan Weiss]

As I asked you at the beginning, did you notice that f(x) is always 0? What makes you think that f(x) is ever nonzero? Did you try to plot it?

x = 0:100;
myf = zeros(size(x));
for F = 0:100
    myf(F) = parameterfun(F);
end
plot(x,myf)

Alan Weiss

MATLAB mathematical toolbox documentation

[Muna Shehan]

Thank Mr. Alan Weiss for your reply. I noticed the objective function is taken (0) value so I think i need to handle this problem and work out to solve it. Thanks again for your replies.

Answer 2

[kamilya MIMOUNI]

Error in Test4 (line 93)

Control_optimal = fmincon(problem.objective)

function Test4

clear;close all; clc

% For starting

fprintf ( 1, '\n' );

fprintf ( 1, 'Test1:\n' );

fprintf ( 1, ' MATLAB version:\n' );

fprintf ( 1, ' A program to demonstrate the finite element method.\n' );

% Geometry data

rin = 0.2; rex = 1;

gd1 = [1;0;0;rex];

gd2 = [1;0;0;rin];

gd = [gd1,gd2];

ns = (char('Cext','Cint'))';

sf = 'Cext - Cint';

[dl, bt] = decsg(gd, sf, ns);

model = createpde;

geometryFromEdges(model,dl);

pdegplot(model,'EdgeLabels','on','FaceLabels','on')

msh = generateMesh(model,'Hmax',0.05,'GeometricOrder','linear','Jiggle','on');

[p,e,t] = meshToPet(msh);

% [p,e,t] = refinemesh(dl,p,e,t);

coordinates = p';

elements3 = t(1:3,:)';

Ne = findNodes(msh,'region','Edge',(1:4));

Ni = findNodes(msh,'region','Edge',(5:8));

neumann = [Ni(:),[Ni(2:end)';Ni(1)]];

dirichlet = [Ne(:),[Ne(2:end)';Ne(1)]];

BoundNodes = unique(dirichlet);

% Export Initial Data

load flux_data.mat flux_data resultss

uintrp = interpolateSolution(resultss,p(1,:),p(2,:));

uintrp(Ne) = 1;

[uintrpx,uintrpy] = pdegrad(p,t,uintrp); % au centre des triangles

uintrpxn = pdeprtni(p,t,uintrpx); % au noeuds

uintrpyn = pdeprtni(p,t,uintrpy); % au noeuds

Nx = @(z)2*z(1,:);

Ny = @(z)2*z(2,:);

z = p(:,Ne);

nx = Nx(z); nx = nx(:);

ny = Ny(z); ny = ny(:);

epsilon = 3/10;

h_flux = 1./(sqrt(nx.^2+ny.^2).*(1+epsilon*z(1,:)')).*(uintrpxn(Ne).*nx+uintrpyn(Ne).*ny);

% Direct problem

RHS_step = zeros(size(elements3,1));

a_x = zeros(size(elements3,1));

c_x = zeros(size(elements3,1));

Control = zeros(size(neumann,1),1);

for j=1:size(elements3,1)

xc = sum(coordinates(elements3(j,:),:))/3;

RHS_step(j) = 0;

c_x(j) = 1./(1+epsilon*xc(1));

a_x(j) = 0;

end

psi_d = ones(size(BoundNodes,1),1); % fct g sur Gamma_ex

for j=1:size(neumann,1)

x_m = sum(coordinates(neumann(j,:),:))/2; % vecteur ligne: coordonnées du milieu du segment

Control(j) = 1; % le controle v sur Gamma_in

end

[psi,psi_flux] = Direct_Problem(elements3,neumann,coordinates,c_x,a_x,RHS_step,Control,psi_d,BoundNodes,p,t, ...

epsilon,nx,ny,Ne,z);

% return

% Cost function

regu = 1/1000;

Integrand = psi_flux - h_flux;

Cost_J = Cost_Function(Control,Integrand,coordinates,Ne,Ni,rex,rin,regu);

if(Cost_J==0)

disp('Your intial guess v = 1 is exactely your solution')

disp('There is no need for optimization')

end

% disp('fmincon start')

% options = optimoptions('fmincon','Display','iter','Algorithm','sqp','PlotFcns', ...

% {@optimplotfval,@optimplotfirstorderopt},'TolFun',4.0e-4);

% options = optimoptions(@fmincon,'GradObj','on','Display','iter','TolFun',1.0e-4, ...

% 'PlotFcns',{@optimplotfval,@optimplotfirstorderopt});

lb = 4.0e-2*ones(size(Control,1),1);

ub = 5*ones(size(Control,1),1);

% return

options = optimoptions('fmincon','Display','iter','Algorithm','sqp');

% return

problem.options = options;

% return

problem.solver = 'fmincon';

% return

problem.objective = @(u_epsilon)Moncef_Cost_Function(u_epsilon,elements3,neumann,coordinates,c_x,a_x,RHS_step, ...

BoundNodes,p,t,epsilon,nx,ny,Ne,Ni,z,rex,rin,regu,h_flux)

% return

problem.x0 = Control;

problem.nonlcon = [];

Control_optimal = fmincon(problem.objective)

return

Control_optimal = fmincon(@(u_epsilon)Moncef_Cost_Function(u_epsilon,elements3,neumann,coordinates,c_x,a_x,RHS_step, ...

BoundNodes,p,t,epsilon,nx,ny,Ne,Ni,z,rex,rin,regu,h_flux), ...

Control,[],[],[],[],lb,ub,[],options);

% Retour au probleme direct final

[psi,psi_flux] = Direct_Problem(elements3,neumann,coordinates,c_x,a_x,RHS_step,Control_optimal,psi_d,BoundNodes,p,t, ...

epsilon,nx,ny,Ne,z);

ep = 0.2;

Nb = findNodes(msh,'box',[-1-ep -1],[-0.5 0.5]);

%[i,c] = pdesdp(p,e,t,sdl);

figure(10)

pdemesh(model_initial)

hold on

plot(msh.Nodes(1,Nb),msh.Nodes(2,Nb),'or','MarkerFaceColor','g')

psi_full = full(psi);

psi_max = max(psi_full(Nb));

indice2 = find(psi_full<=psi_max+0.07 & psi_full>=psi_max-0.07);

k = convhull(p(1,indice2),p(2,indice2));

figure(12)

pdemesh(model)

hold on

plot(p(1,k),p(2,k),'r-')

figure(11)

pdegplot(model)

hold on

plot(msh.Nodes(1,indice2),msh.Nodes(2,indice2),'.r','MarkerFaceColor','r')

% return

[i,c] = pdesdp(p,e,t,sdl) % given mesh data p, e, and t and a list of subdomain numbers sdl, the function returns all points belonging to those subdomains. A point can belong to several subdomains, and the points belonging to the domains in sdl are divided into two disjoint sets. i contains indices of the points that wholly belong to the subdomains listed in sdl, and c lists points that also belongs to the other subdomains.

c = pdesdp(p,e,t,sdl) % returns indices of points that belong to more than one of the subdomains in sdl.

i = pdesdt(t,sdl) % given triangle data t and a list of subdomain numbers sdl, i contains indices of the triangles inside that set of subdomains.

i = pdesde(e,sdl) % given edge data e, it extracts indices of outer boundary edges of the set of subdomains.

%====================================

[p,~,t] = meshToPet(mesh);

%Get the right area for integration for all named areas in sdl

tsub = t(:,t(4,:)==sdl(1));

%Get the interpolated midpoint

ut = pdeintrp(p,tsub,u);

[a,~,~,~] = pdetrg(p,tsub);

a_sum=sum(a);

%Get mean value of solution in the area

u_mean = sum(ut.*a)/a_sum;

%===========================================

Answer 3

[kamilya MIMOUNI]

please help me guys