Optimization problem - GlobalSearch is not generating initial vector in MATLAB

Question

INTRO

I have to solve optimization problem, containing a set of ordinary differential equations. Whats more I need to find global minimum. So far I've only managed to find few local minimums, but examining manually which one is global is tricky for so complex problem.

Description

Let's say I have 3 Matlab files:

1. Main.m - its the general script file, which is a driver for all others
2. ODE_set.m - this hold all my ODE's system
3. Opt_query.m - this contain the definition of optimization objectives, and all equations needed to determined the minimum. Also from this file the ODE solver is started.

I can look for local minimum with this code, and it works totally fine, without any errors:

Main.m

global k1 k2 k3
x0 = [0.5 8 13];
lu =[ 0.9 25.74 60.9];
lb =[ 0.1 0.74 0.9];
[k,fval,exitflag,output] = fmincon('Opt_query',x0,[],[],[],[],lb,lu);
disp(k1)
disp(k2)
disp(k3)

However as I mentioned above I need global minimum so I started to try with GlobalSearch. Here is my code:

x0 = [0.5 8 13];
lu =[ 0.9 25.74 60.9];
lb =[ 0.1 0.74 0.9];
gs = GlobalSearch;
problem = createOptimProblem('fmincon',Opt_query,x0,[],[],[],[],lb,lu);
[x,fval,exitflag,output] = run(gs,problem);

It ends with an error:

Not enough input arguments.

Error in Opt_query (line 3)

k1 = p(1)

Here is the structure of Opt_query.m:

function q = Opt_query(p)
    global k1 k2 k3
    k1 = p(1)
    k2 = p(2)
    k3 = p(3)
    y0=[1 2 3]
    [t,y] = ode45(@ODE_set,[0, 130],y0,k1,k2,k3)
    t_steps=0:5:130;
    y_steps=interp1(t,y,t_steps);
    for j = 1:27
        Z = Z+ abs(y_steps(j,1)-5)+abs(y_steps(j,2)-7)+abs(y_steps(j,3)-9);
    end
    q = sqrt(Z);
end

Here is the structure of ODE_set.m:

function dydt = ODE_set(t, y, k1, k2, k3)

F = 20.1;
A_in = 2.5;
V = 100;
k = 0.150; 
A = y(1);
B = y(2);
C = y(3);

n = numel(y);
dydt = zeros(n,1);
dydt(1) = F/V*(A_in - A) - k1*A^2;
dydt(1) = F/V*(B) - k2*B^2;
dydt(1) = F/V*(C) - k3*C^2;

The Problem

So my question is - why when I use fmincon alone, it works fine, but when initiated with GlobalSearch it doesn't generate initial p vector which holds parameters to optimized? Did I miss something, or maybe GlobalSearch require different syntax?

Answer 1

Few mistakes

---

In Opt_query change this

[t,y] = ode45(@ODE_set,[0, 130],y0,k1,k2,k3);

to this

[t,y] = ode45(@(t,y)ODE_set(t,y, k1, k2, k3),[0, 130],y0);

After tspan = [0,30] you should only use the initial conditions, do not include k1, k2, k3

• Just take the function ODE_set which has 5 inputs t,y, k1, k2, k3 make a new one with two inputs t , ythen k1, k2, k3 will be kept as fixed inputs

new_ODE_set = @(t, y)ODE_set(t, y, k1, k2, k3)

• Here is the syntax for ode45

[t,y] = ode45(odefun,tspan,y0)

odefun is a function of t, y only

• Since new_ODE_set is a function of t, y as well it can be used now as odefun

• but actually we still have k1, k2, k3 inside it

---

for j = 1:27
   Z = Z+ abs(y_steps(j,1)-5)+abs(y_steps(j,2)-7)+abs(y_steps(j,3)-9);
end

Here you should initialize Z to zero before running the loop, like this

Z = 0;
for j = 1:27
    Z = Z+ abs(y_steps(j,1)-5)+abs(y_steps(j,2)-7)+abs(y_steps(j,3)-9);
end

---

In ODE_set you're only using index 1 for dydt

Change this

dydt(1) = F/V*(A_in - A) - k1*A^2;
dydt(1) = F/V*(B) - k2*B^2;
dydt(1) = F/V*(C) - k3*C^2;

to this

dydt(1) = F/V*(A_in - A) - k1*A^2;
dydt(2) = F/V*(B) - k2*B^2;
dydt(3) = F/V*(C) - k3*C^2;

---

To sum up

---

• Main.m

global k1 k2 k3
x0 = [0.5 8 13];
lu =[ 0.9 25.74 60.9];
lb =[ 0.1 0.74 0.9];
gs = GlobalSearch;

problem = createOptimProblem('fmincon','objective',...
    @Opt_query,'x0',x0,'lb',lb,'ub',lu);

[x,fval,exitflag,output] = run(gs,problem);
disp(k1)
disp(k2)
disp(k3)

---

• ODE_set.m

function dydt = ODE_set(t, y, k1, k2, k3)

    F = 20.1;
    A_in = 2.5;
    V = 100;
    k = 0.150; 
    A = y(1);
    B = y(2);
    C = y(3);

    n = numel(y);
    dydt = zeros(n,1);
    dydt(1) = F/V*(A_in - A) - k1*A^2;
    dydt(2) = F/V*(B) - k2*B^2;
    dydt(3) = F/V*(C) - k3*C^2;
end

---

• Opt_query.m

function q = Opt_query(p)


    global k1 k2 k3
    k1 = p(1);
    k2 = p(2);
    k3 = p(3);
    y0=[1 2 3];



    [t,y] = ode45(@(t,y)ODE_set(t,y, k1, k2, k3),[0, 130],y0);
    t_steps=0:5:130;
    y_steps=interp1(t,y,t_steps);
    Z = 0;
    for j = 1:27
        Z = Z+ abs(y_steps(j,1)-5)+abs(y_steps(j,2)-7)+abs(y_steps(j,3)-9);
    end
    q = sqrt(Z);
end