Using global optimization Matlab toolbox and genetic algorithm to minimize black-box function

I have following function:

function w = objective_function(a,b,c,d)
    function w = rh(t,y) 
        w = [y(2);6*cos(9*t) - 10*sin(3*t) - a*y(1) - b*(1-y(1)^2)*y(2)];
    end
    t = linspace(0,20,20000);
    [t,y] = ode45(@rh,t,[c;d]);
    u = 2*sin(3*t);
    w = norm(y(:,1)-u,Inf)/norm(u,Inf);
end

and I would like to find its minimum using genetic algorithm and globaloptimization toolbox, so I wrote following script:

a = optimvar("a","LowerBound",0,"UpperBound",10);
    b = optimvar("b","LowerBound",0,"UpperBound",10);
    c = optimvar("c","LowerBound",-5,"UpperBound",5);
    d = optimvar("d","LowerBound",-5,"UpperBound",5);
    prob = optimproblem("Objective",objective_function(a,b,c,d));
    options = optimoptions("ga","PlotFcn","gaplotbestf");
    rng default
    [sol,fval] = solve(prob,"Solver","ga","Options",options)

and I got following error when I run it:

     Error using superiorfloat
    Inputs must be floats, namely single or double.

    Error in odearguments (line 114)
        dataType = superiorfloat(t0,y0,f0);

    Error in ode45 (line 104)
        odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, 
        tspan, y0, options, varargin);

    Error in objective_function (line 13)
        [t,y] = ode45(rh,[0,10],[c;d]);

    Error in script (line 6)
        prob = optimproblem("Objective",objective_function(a,b,c,d));

How to fix this problem?

Answer 1

[Torsten]

rng default
a = optimvar("a","LowerBound",0,"UpperBound",10);
b = optimvar("b","LowerBound",0,"UpperBound",10);
c = optimvar("c","LowerBound",-5,"UpperBound",5);
d = optimvar("d","LowerBound",-5,"UpperBound",5);
obj = fcn2optimexpr(@objective_function,a,b,c,d)
prob = optimproblem("Objective",obj);
%[x0.a,x0.b,x0.c,x0.d] = deal(1,1,1,1);
options = optimoptions("ga","PlotFcn","gaplotbestf");
show(prob)
[sol,fval] = solve(prob,"Solver","ga","Options",options)
%[sol,fval] = solve(prob,'x0',x0)
    
function w = objective_function(a,b,c,d)
    function w = rh(t,y) 
        w = [y(2);6*cos(9*t) - 10*sin(3*t) - a*y(1) - b*(1-y(1)^2)*y(2)];
    end
    t = linspace(0,20,20000);
    [t,y] = ode45(@rh,t,[c;d]);
    u = 2*sin(3*t);
    w = norm(y(:,1)-u,Inf)/norm(u,Inf);
end 

[Pawel Labedzki]

Thank you!