fmincon converges to initial value

Context:

The problem I'm solving is in finding the optimal placement and size of a piezoelectric patch on a beam such that the modal force will be maximized. In this code I calculate the modal shapes using the Ritx method, and then apply an equation to get the modal force and then sum over the different modes and divide by the number of modes I'm adding. The goal of the constraint is to minimize the drifting of the natural frequencies (i.e. imagine a very thin and short beam with a large piezoeelctric patch - not good). The way I implement this is by constricting all modified natural frequencies to be within a n-dimensional ball of radius c*magnitude(omega), where c is a number from 0 to 1.

Is this smooth?

From the below plots, it looks like the function seems to be quite smooth. The left plot is a color map showing the constraint function, and the right plot is a color map of the function we want to minimize.

Optimization Code:

% Point Density
M = 1;

H = linspace(0.25e-3,0.24e-2,M);
E = 200*10^9;
rho = 7700;

xi_opt = zeros(2,M);

for j = 1:M
    [xi_opt(:,j)] = Optimizer(H(j),E,rho);
end

function [xi1] = Optimizer2(h,E,rho)



end

function [xi] = Optimizer(h,E,rho)

% Constraint value
c = 0.1;
L = 1;

xi_0 = [0.2;0.34];
A = [1,2];
b = [L];
Aeq = [];
beq = [];
lb = eps.*ones(2,1);
ub = Inf.*ones(2,1);
options = optimoptions('fmincon','display','iter');

xi = fmincon(@(xi) f(xi,h,E,rho),xi_0,A,b,Aeq,beq,lb,ub,@(xi) ...
    nonlcon(xi,c,h,E,rho),options);
end

function [ciq,ceq] = nonlcon(xi,c,h,E,rho)
[~,g] = f(xi,h,E,rho);
ciq = g-c;
ceq = [];
end

function [metric1,metric2] = f(xi,h,E,rho)
%% Material Properties

% Beam properties
L=1;        % beam length (m)
b=0.1;	    % beam width (m)

%{
h=.24e-3;   % thickness (m)
E=200*10^9; % Young's modulus
rho=7700;   % density
%}

% Properties of the Piezoelectric Film
E_p=30.336e9;   %Young's modulus (Pa)
d33=460e-12;    %piezo constant, elongating type (C/N)
h_p=0.30e-3;    %thickness of MFC, (m)
l_p=0.45e-3;    %distance between poles (m)
rho_p=5440;     %density (kg/m^3)

%% Modal settings
N=7;      % number of Ritz vectors to use (with 13, first five non-zero 
%          % modes are well converged)
epw=100;   % number of elements per half-wave to discretize beam into

% Voltage amplitude
V=400;          %Amplitude (+/-400 is about the maximum safe limit)

%% Piezoelectric Calcs

% Distance from neutral axis to geometric mid-plane 
y0=(E_p*h_p*(h+h_p))/(2*(E*h+E_p*h_p));
% Electromechanical coupling coefficient
chi=(d33*E_p*b*h_p/l_p)*((.5*(h+h_p)-y0));

% No adjustments for patch stiffness or mass
[wn, ~, ~]=DDSL_Beam_Modes(L,b,h,E,rho,N,epw);

% Adjustments for patch stiffness or mass
[wnp, Phip, nodemap]=DDSL_Beam_Modes_p(L,b,h,h_p,E,E_p,rho,...
    rho_p,N,epw,xi);

%% Calcs on modal outputs

% Initialize secondary counters for average and distance
g = 0;
h = 0;
metric1 = 0;
metric2 = 0;

% Initialize modal force vector
Q = zeros(N,1);

% differential distance
dx = nodemap(2)-nodemap(1);

% Node points at beginning and end of piezoelectric element
k = dsearchn(nodemap',[xi(2)*L; L*(xi(1)+xi(2))]);

% Compute modal force (See Erturk and Inman Book Eq. 3.31)
for n = 1:N
    dphi=gradient(Phip(:,n))/dx;         % Use piezo-adjusted modes
    Q(n)=chi*V*(dphi(k(2))-dphi(k(1)));

    if (n > 2)&&(n < 8)
        metric1 = metric1+abs(Q(n));
        metric2 = metric2+((wn(n)-wnp(n)))^2;
        g = g+1;
        h = h+(wn(n))^2;
    end
end

% Normalize metrics
metric1 = -metric1/g;
metric2 = sqrt((metric2/h));

end

Upon request, I can provide code for the script that produced the graph as well as the two functions that this code internally references, but I honestly don't think that those fuctions are the source of the problem here.

Question:

When I use different initial conditions, I notice that fmincon seems to get stuck in a local min very near the initial conditions - one that clearly isn't correct.

When I set M = 1, and use the initial condition xi = [0.2;0.34], I get

xi_opt

xi_opt =

   0.197630581341251
   0.340996028109980

>> 

I also used some options that would display the function value at each step. The following is the command window output:

DDSL_Piezo_Optimizer
                                            First-order      Norm of
 Iter F-count            f(x)  Feasibility   optimality         step
    0       3   -4.373705e+00    0.000e+00    8.139e-02
    1      15   -4.373943e+00    0.000e+00    8.420e-02    1.700e-03
    2      26   -4.374048e+00    0.000e+00    8.543e-02    7.381e-04
    3      39   -4.374060e+00    0.000e+00    8.558e-02    8.046e-05
    4      50   -4.374065e+00    0.000e+00    8.563e-02    3.519e-05
    5      61   -4.374067e+00    0.000e+00    8.564e-02    1.539e-05
    6      76   -4.374067e+00    0.000e+00    8.564e-02    4.208e-07
    7      87   -4.374067e+00    0.000e+00    8.564e-02    1.841e-07
    8     100   -4.374067e+00    0.000e+00    8.567e-02    2.014e-08
    9     110   -4.374067e+00    0.000e+00    8.561e-02    1.762e-08
   10     124   -4.374067e+00    0.000e+00    8.565e-02    9.636e-10
   11     136   -4.374067e+00    0.000e+00    8.563e-02    2.108e-10
   12     146   -4.374067e+00    0.000e+00    8.564e-02    1.844e-10

Note that f(x) converged to a value of -4.37, yet just from the color bar on Figure 1, we know that the minimum must be near -14.

Although it seems like this function is smooth it appears like it has a very large number of local minima. Have I made some mistakes causing this? If fmincon isn't the correct algorithm to be using, then what would you reccomend? What trouble shooting guidelines would you give to help with solving these issues.

Answer 1

[Alan Weiss]

I haven't looked at your problem in detail, but I suppose that the issue is related to those found in Optimizing a Simulation or ODE.

Good luck,

Alan Weiss

MATLAB mathematical toolbox documentation