MATLAB Answers

Obtained results are far from the specified objective function when using fmincon with a Multi-Start.

51 views (last 30 days)
Hasan Hassoun
Hasan Hassoun on 23 Jan 2021
Commented: Hasan Hassoun on 23 Jan 2021
This question was flagged by John D'Errico
The obtained results are far from the specified objective function when using fmincon with Multi start.
Any comment/help is welcome, but please read the question and code well before any answer!
function ObjectiveFunction = N_var(x0)
N=5;
i=1:1:N; o=N+1:1:2*N; p=((2*N)+1):1:3*N;
Ln(i) = x0(i); % Length of the neck
S(i) = x0(o); % Area of the neck
V(i) = x0(p); % Volume of cavity
c0=340.3; % Speed of sound in the air
density=1.225; % Air density
syms w;
% Variables
Lv=2.5e-3; %length of the cicular cavity
Rn(i)= (S(i)./pi()).^0.5; % Radius of the circular neck
Rv(i)=(V(i)./(Lv.*pi())).^0.5; % Radius of the circular cavity
% Calculations
% Effective neck length
q=2:1:N;g=2:1:N-1;
Leff(1)=Ln(1)+(0.82.*(Rn(1)./2))+(0.85.*Rn(1).*(1-(1.25.*(Rn(1)./Rv(1)))));
Leff(q)=Ln(q)+(0.85.*Rn(q).*(1-(1.25.*(Rn(q)./Rv(q-1)))))+(0.85.*Rn(q).*(1-(1.25.*(Rn(q)./Rv(q)))));
Mass=[];
Mass(1)=density.*S(1).*Leff(1);
Mass(q)=density.*S(q).*Leff(q);
% Matrix input
M=zeros(N,N);
M =subs(M,M(1,1),((-1.*w.^2.*Mass(1))+((density.*(c0.^2).*((S(1)).^2))./V(1))));
M(1,2)= (-1.*density.*(c0.^2).*((S(1)).^2).*((S(2)).^2))./V(1);
M(2,1)= (-1.*density.*(c0.^2).*((S(1)).^2).*((S(2)).^2))./V(1);
for q=2:1:N
M(q,q)=((-1.*w.^2.*Mass(q))+((density.*(c0.^2).*((S(q)).^2)).*((1./V(q-1))+(1./V(q)))));
end
for g=2:1:N-1
M(g+1,g)= (-1.*density.*(c0.^2).*((S(g)).^2).*((S(g+1)).^2))./V(g);
M(g,g+1)= (-1.*density.*(c0.^2).*((S(g)).^2).*((S(g+1)).^2))./V(g);
end
M=M.*tril( triu(true(size(M)),-1), +1); % create the zeros of the matrix at upper and lower limits(diag+1 until n and diag-1 until n) %% 2nd option (% for k=1:1:N-2 % M(k,(k+2):N)=0; % M((k+2):N,k)=0; % end)
detM=det(M);
% Results
w = vpasolve(detM,w);
%Because there are no parameters in this solution, use vpa to approximate it numerically.
F_resonance =double(w./(2*pi()));
F_resonance=sort(F_resonance);
speed=3000/60; % fundamental frequency (rpm/60);
F_Harmonic(i)= (i/2).*speed;
F_Harmonic= (F_Harmonic)'; % Transpose F_Harmonic
ObjectiveFunction = sum(abs (F_resonance(o)-F_Harmonic(i)));
end
%% Optimization code
clear all; clc;
N=5; %Number of DOF;
%% Optimization function
A=[]; b=[]; Aeq=[]; beq=[]; nonlcon=[]; %b and beq are vectors,A and Aeq are matrices,nonlcon are the nonlinear constraints.Set arguments to [] when there are no constraints.
lb=repmat(1e-5,1,3*N); ub=repmat(2,1,3*N); % lb and ub are the vectors of lower and upper bounds of the solution.
x0= lb + rand(size(lb)).*(ub - lb); % is the initial point for solution % x0(1) = Ln1; % x0(2) = Ln2; % x0(3) = S1; % x0(4) = S2; % x0(5) = V1; % x0(6) = V2;
figure (1) % plotting
options = optimset('MaxFunEvals',Inf,'MaxIter',10000,'Display','iter','PlotFcn', {@optimplotfval});
[solution,objectivefunction,exitflag,output] = fmincon(@N_var,x0,A,b,Aeq,beq,lb,ub,nonlcon,options); %"fval" objective function value at solution,"exitflag" reason fmincon stopped.
%% TL without Multi-start(ms)/Global search(gs)
i=1:1:N; o=N+1:1:2*N; p=((2*N)+1):1:3*N;
Ln(i) = solution(i); % Length of the neck
S(i) = solution(o); % Area of the neck
V(i) = solution(p); % Volume of cavity
c0=340.3; % Speed of sound in the air
density=1.225; % Air density
Sd=4.3*4.3e-4; % duct area
Lv=2.5e-3; %length of the cicular cavity
freq=0:0.25:500; omega=2.*pi().*freq;
speed=3000/60; % fundamental frequency (rpm/60);
F_Harmonic(i)= (i/2).*speed;
F_Harmonic= (F_Harmonic)'; % Transpose F_Harmonic
syms w;
Rn(i)= (S(i)./pi()).^0.5; % Radius of the circular neck
Rv(i)=(V(i)./(Lv.*pi())).^0.5; % Radius of the circular cavity
% % Calculations
% Effective neck length
q=2:1:N;g=2:1:N-1;
Leff(1)=Ln(1)+(0.82.*(Rn(1)./2))+(0.85.*Rn(1).*(1-(1.25.*(Rn(1)./Rv(1)))));
Leff(q)=Ln(q)+(0.85.*Rn(q).*(1-(1.25.*(Rn(q)./Rv(q-1)))))+(0.85.*Rn(q).*(1-(1.25.*(Rn(q)./Rv(q)))));
Mass=[];
Mass(1)=density.*S(1).*Leff(1);
Mass(q)=density.*S(q).*Leff(q);
% Matrix input
M=zeros(N,N);
M =subs(M,M(1,1),((-1.*w.^2.*Mass(1))+((density.*(c0.^2).*(S(1).^2))./V(1))));%[Mass(1),density,c0,S(1),V(1)],[Mass(1),density,c0,S(1),V(1)]);
M(1,2)= (-1.*density.*(c0.^2).*(S(1).^2).*(S(2).^2))./V(1);
M(2,1)= (-1.*density.*(c0.^2).*(S(1).^2).*(S(2).^2))./V(1);
%
for q=2:1:N
M(q,q)=((-1.*w.^2.*Mass(q))+((density.*(c0.^2).*(S(q).^2)).*((1./V(q-1))+(1./V(q)))));
end
for g=2:1:N-1
M(g+1,g)= (-1.*density.*(c0.^2).*(S(g).^2).*(S(g+1).^2))./V(g);
M(g,g+1)= (-1.*density.*(c0.^2).*(S(g).^2).*(S(g+1).^2))./V(g);
end
%
M=M.*tril( triu(true(size(M)),-1), +1); % create the zeros of the matrix at upper and lower limits(diag+1 until n and diag-1 until n) %% 2nd option (% for k=1:1:N-2 % M(k,(k+2):N)=0; % M((k+2):N,k)=0; % end)
inverse_M =inv(M);
% %Results
Z=(1./(1i.*w.*density.*c0.*S(1).*inverse_M(1,1)));
TL=20*log10(abs(1+(S(1)./(2*Sd.*Z))));
TL=subs(TL,w,omega);
TL=double(TL);
%% TL based on Multi-start(ms)
figure(1)
opts_ms = optimoptions(@fmincon);
problem_ms = createOptimProblem('fmincon','objective',@N_var,'x0',x0,'lb',lb,'ub',ub,'options',opts_ms);
[solution_ms,objectivefunction_ms,exitflag_ms,iteration_ms,output_ms] = run(MultiStart('PlotFcn',@gsplotbestf),problem_ms,100)
Ln_ms(i) = solution_ms(i); % Length of the neck
S_ms(i) = solution_ms(o); % Area of the neck
V_ms(i) = solution_ms(p); % Volume of cavity
Rn_ms(i)= (S_ms(i)./pi()).^0.5; % Radius of the circular neck
Rv_ms(i)=(V_ms(i)./(Lv.*pi())).^0.5; % Radius of the circular cavity
% Calculations
% Effective neck length
Leff_ms(1)=Ln_ms(1)+(0.82.*(Rn_ms(1)./2))+(0.85.*Rn_ms(1).*(1-(1.25.*(Rn_ms(1)./Rv_ms(1)))));
Leff_ms(q)=Ln_ms(q)+(0.85.*Rn_ms(q).*(1-(1.25.*(Rn_ms(q)./Rv_ms(q-1)))))+(0.85.*Rn_ms(q).*(1-(1.25.*(Rn_ms(q)./Rv_ms(q)))));
Mass_ms=[];
Mass_ms(1)=density.*S_ms(1).*Leff_ms(1);
Mass_ms(q)=density.*S_ms(q).*Leff_ms(q);
% Matrix input
M_ms=zeros(N,N);
M_ms =subs(M_ms,M_ms(1,1),((-1.*w.^2.*Mass_ms(1))+((density.*(c0.^2).*(S_ms(1).^2))./V_ms(1))));%[Mass(1),density,c0,S(1),V(1)],[Mass(1),density,c0,S(1),V(1)]);
M_ms(1,2)= (-1.*density.*(c0.^2).*(S_ms(1).^2).*(S_ms(2).^2))./V_ms(1);
M_ms(2,1)= (-1.*density.*(c0.^2).*(S_ms(1).^2).*(S_ms(2).^2))./V_ms(1);
for q=2:1:N
M_ms(q,q)=((-1.*w.^2.*Mass_ms(q))+((density.*(c0.^2).*(S_ms(q).^2)).*((1./V_ms(q-1))+(1./V_ms(q)))));
end
for g=2:1:N-1
M_ms(g+1,g)= (-1.*density.*(c0.^2).*(S_ms(g).^2).*(S_ms(g+1).^2))./V_ms(g);
M_ms(g,g+1)= (-1.*density.*(c0.^2).*(S_ms(g).^2).*(S_ms(g+1).^2))./V_ms(g);
end
M_ms=M_ms.*tril( triu(true(size(M_ms)),-1), +1); % create the zeros of the matrix at upper and lower limits(diag+1 until n and diag-1 until n) %% 2nd option (% for k=1:1:N-2 % M(k,(k+2):N)=0; % M((k+2):N,k)=0; % end)
inverse_M_ms =inv(M_ms);
%Results
Z_ms=(1./(1i.*w.*density.*c0.*S_ms(1).*inverse_M_ms(1,1)));
TL_ms=20*log10(abs(1+(S_ms(1)./(2*Sd.*Z_ms))));
TL_ms=subs(TL_ms,w,omega);
TL_ms=double(TL_ms);
%% TL based on Global search(gs)
opts_gs = optimoptions(@fmincon);
problem_gs = createOptimProblem('fmincon','objective',@N_var,'x0',x0,'lb',lb,'ub',ub,'options',opts_gs);
gs=GlobalSearch (MultiStart('PlotFcn',@gsplotbestf));% or use (MultiStart('PlotFcn',@gsplotbestf)) instead of ('Display','iter')
[solution_gs,objectivefunction_gs] = run(gs,problem_gs)
Ln_gs(i) = solution_gs(i); % Length of the neck
S_gs(i) = solution_gs(o); % Area of the neck
V_gs(i) = solution_gs(p); % Volume of cavity
Rn_gs(i)= (S_gs(i)./pi()).^0.5; % Radius of the circular neck
Rv_gs(i)=(V_gs(i)./(Lv.*pi())).^0.5; % Radius of the circular cavity
% % Calculations
% Effective neck length
Leff_gs(1)=Ln_gs(1)+(0.82.*(Rn_gs(1)./2))+(0.85.*Rn_gs(1).*(1-(1.25.*(Rn_gs(1)./Rv_gs(1)))));
Leff_gs(q)=Ln_gs(q)+(0.85.*Rn_gs(q).*(1-(1.25.*(Rn_gs(q)./Rv_gs(q-1)))))+(0.85.*Rn_gs(q).*(1-(1.25.*(Rn_gs(q)./Rv_gs(q)))));
Mass_gs=[];
Mass_gs(1)=density.*S_gs(1).*Leff_gs(1);
Mass_gs(q)=density.*S_gs(q).*Leff_gs(q);
% % Matrix input
M_gs=zeros(N,N);
M_gs =subs(M_gs,M_gs(1,1),((-1.*w.^2.*Mass_gs(1))+((density.*(c0.^2).*(S_gs(1).^2))./V_gs(1))));%[Mass(1),density,c0,S(1),V(1)],[Mass(1),density,c0,S(1),V(1)]);
M_gs(1,2)= (-1.*density.*(c0.^2).*(S_gs(1).^2).*(S_gs(2).^2))./V_gs(1);
M_gs(2,1)= (-1.*density.*(c0.^2).*(S_gs(1).^2).*(S_gs(2).^2))./V_gs(1);
for q=2:1:N
M_gs(q,q)=((-1.*w.^2.*Mass_gs(q))+((density.*(c0.^2).*(S_gs(q).^2)).*((1./V_gs(q-1))+(1./V_gs(q)))));
end
for g=2:1:N-1
M_gs(g+1,g)= (-1.*density.*(c0.^2).*(S_gs(g).^2).*(S_gs(g+1).^2))./V_gs(g);
M_gs(g,g+1)= (-1.*density.*(c0.^2).*(S_gs(g).^2).*(S_gs(g+1).^2))./V_gs(g);
end
M_gs=M_gs.*tril( triu(true(size(M_gs)),-1), +1); % create the zeros of the matrix at upper and lower limits(diag+1 until n and diag-1 until n) %% 2nd option (% for k=1:1:N-2 % M(k,(k+2):N)=0; % M((k+2):N,k)=0; % end)
inverse_M_gs =inv(M_gs);
%Results
Z_gs=(1./(1i.*w.*density.*c0.*S_gs(1).*inverse_M_gs(1,1)));
TL_gs=20*log10(abs(1+(S_gs(1)./(2*Sd.*Z_gs))));
TL_gs=subs(TL_gs,w,omega);
TL_gs=double(TL_gs);
% TL plotting
figure (1)
plot(freq,TL,'red',freq,TL_ms,'black',freq,TL_gs,'blue')
xlabel('Frequency (Hz)'); ylabel('Transmission Loss (dB)');
legend('without Multi-start','Multi start','Global search')
solution, objectivefunction
solution_ms, objectivefunction_ms
solution_gs, objectivefunction_gs
  4 Comments
Hasan Hassoun
Hasan Hassoun on 23 Jan 2021
I posted the code again. When using fmincon without multi-start, the objective function converged, and the obtained results fit well for what is specified in the objective function. However, when using multi-start, the obtained results are far from what I expect (Please check the attached figure above).

Sign in to comment.

Answers (1)

John D'Errico
John D'Errico on 23 Jan 2021
Edited: John D'Errico on 23 Jan 2021
No numerical optimization that applies to a completely general nonlinear function will converge to the global soution regardless of the starting points.
A global search methods like multi-start uses randomization with multiple start points, thus just like it sounds. But that means it starts the optimizer out at randomly chosen different start points, hoping that one of them will land in the basin of attraction of the global optimum. However, since the start points are randomly chosen, and the basin of attraction of the global optimizer may be infinitely small, there is no assurance they will always yield the same result.
If you want a better success rate, then you need to use a deeper more intensive search, so more start points. But you will never have a guaranteed success rate of finding the same global optimum if your function is nasty enough.
Similar conclusions apply to any such scheme that uses randomized methods.
And any scheme that does not use randomized methods to overcome the limitations of bad starting values will suffer the fate of poor starting values. You end up with a non-global solution. But it will be repeatable. Is that "better"? Is a "reliable" poor solution (that depends completely on the qualty of your starting values) better than a sometimes unreliable solution? If it did not converge to the global optimum, but you do not know it did not, would you honestly call that "reliable"? Are one of those people who thinks if a tree falls in the woods and nobody is there to hear it, there was no sound?
Best is to understand the tools you are using. Here, I would suggest you learn about things like basins of attraction as they apply to numerical optimizers.

Community Treasure Hunt

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

Start Hunting!