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

조회 수: 7 (최근 30일)
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>

채택된 답변

Alan Weiss
Alan Weiss 2016년 6월 8일
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
  댓글 수: 5
Alan Weiss
Alan Weiss 2016년 6월 10일
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
Muna Shehan 2016년 6월 11일
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.

댓글을 달려면 로그인하십시오.

추가 답변 (2개)

kamilya MIMOUNI
kamilya MIMOUNI 2019년 10월 26일
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;
%===========================================

kamilya MIMOUNI
kamilya MIMOUNI 2019년 10월 26일
please help me guys

카테고리

Help CenterFile Exchange에서 Solver-Based Nonlinear Optimization에 대해 자세히 알아보기

태그

Community Treasure Hunt

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

Start Hunting!

Translated by