Using Matlab's fmincon to find optimal parameter of SDOF system having discrete white noise input

Question

I am trying to use Matlab's fmincon function to find the optimal value of $k$ for the following ODE, which represents a mass-spring system, $$m\ddot{x}(t)+kx(t)=f(t),$$ where $f(t)$ is taken to be constant amplitude white noise. The cost function is taken to be absolute maximum displacement of the mass, that is, $\max_\limits{t}(|x(t)|)$.

The issue I am having is that fmincon takes a very long time to find the optimal value when using a white noise input. When using a sinusoidal input, fmincon finds the optimal value relatively quickly. Note that this is just an illustrative example of a problem I'm experiencing with a more complicated ODE, that is, even with the complicated ODE, having a sinusoidal input the optimal value is found relatively quickly.

First let's generate the white noise input and save it into an .mat file,

stepsize = 0.01;
tgm = 0:stepsize:10;
fgm = randn(numel(tgm),1);
save('gm.mat','tgm','fgm')

The code is broken into three parts:

The objective function and loading the discrete input:

function cost = objfunc(para)
tstart = 0;
tfinal = 10;    % sec

load('gm.mat','tgm','fgm')
options = odeset('RelTol',1e-8 ,'AbsTol',1e-8 );

z0 = [0 ; 0];

[t,z] = ode45(@(t,z) sdof(t,z,para,tgm,fgm), [tstart,tfinal], z0,options);

cost = max(abs(z(:,1)));
end

The ODE:

function dz = sdof(t,z,para,tgm,fgm)
dz = zeros(2,1);
k = para;
m = 1
f = interp1(tgm, fgm, t, 'linear');

dz(1) = z(2);
dz(2) = -(k/m)*z(1) + f/m;
end

Calling fmincon:

para0 = 0.1;
lb = 0;
ub = 1;
para = fmincon(@objfunc, para0, [], [], [], [],lb, ub);

I'm pretty new to using Matlab's Optimization Toolbox so I'm not entirely sure if I'm approaching this problem correctly, or if I should code the white noise differently or input it in a different place.

Thanks for any help or feedback!

Answer 1

Surprisingly, but the problem is not in load. When I run your code with a profiler, it takes 460 seconds, where 396 seconds are spent for interp1, 2808319 calls, that is 141 ms per 1000 calls.

If you open the code of this function, you will find that it makes a lot of things to parse your inputs, then it creates a griddedInterpolant, and then it interpolates. Since we know that we always call it with the same tgm and fgm, let us create the griddedInterpolant on our side, and then pass it to the function as an argument. Ok, now the code with the profiler takes 178 seconds (125 seconds without profiling), where interpolation takes 71 seconds for the same 2808319 calls, that is 25 ms per 1000 calls, more than 5 times faster.

Creating anonymous functions does not change a lot in this case, but I have also updated the code to show how it works.

The function that creates a handle for the cost function:

function obj_fh = get_objfun(tgm, fgm)

    data_set = griddedInterpolant(tgm,fgm);

    tstart = 0;
    tfinal = 10;    % sec
    z0 = [0 ; 0];

    options = odeset('RelTol',1e-8 ,'AbsTol',1e-8 );

    obj_fh = @(para) objfun(para);

        function cost = objfun(para)
            [~,z] = ode45(@(t,z) sdof(t,z,para,data_set), [tstart,tfinal], z0,options);

            cost = max(abs(z(:,1)));
        end
end

The function that implements the ode

function dz = sdof(t,z,para,data_set)
    dz = zeros(2,1);

    k = para;
    m = 1;
    f = data_set(t);

    dz(1) = z(2);
    dz(2) = -(k/m)*z(1) + f/m;
end

The main script

para0 = 0.1;
lb = 0;
ub = 1;

load('gm.mat','tgm','fgm');
objfunc = get_objfun(tgm, fgm);

para = fmincon(objfunc, para0, [], [], [], [],lb, ub);