I use gekko for optimal control of complex network, but if often can't extract hessian matrix (stacked in this step)

Question

I am doing an optimal control problem of complex network of a variation model of SIR model. The network contains 5000 nodes and I will count the number of nodes of each degree (from 4 to 187). The txt file is the file of degrees, the first column is the non-zero degrees in the network and the second line is the number of nodes in this degree.

There are total 76 non-zero degrees, so the number of equations is 76*3. The simulation time is 1000，and my purpose is to find the optimal control. I used IMODE=6 and solver=3.

But the running time is too long. In addition, most of time it can't extract Hessian Matrix, so the running process often stacked. I don't know how to speed up the calculation or make the calculation more stable (slow can be accepted but I need to make sure it is running normally).

Here is the code.

import numpy as np
from gekko import GEKKO
import matplotlib.pyplot as plt
import os

N = 5000
T = 1000
beta = 0.025
phi = 0.005
gamma = 0.01
sigma = 0.001
kmin = 4
kmax = 187
k_mean = 3.9968

data = np.loadtxt("degrees_5000_4_187_72.txt")
pk = data[:,2].tolist()
degrees = data[:,0].tolist()
nums_k = len(degrees)

m = GEKKO(server="https://byu.apmonitor.com")#remote = False

#set the variables
Uk = m.Array(m.MV,nums_k)
for i in range(nums_k):
    Uk[i].LOWER = 0
    Uk[i].UPPER = 1
    Uk[i].STATUS = 1

Sk_initial = 0.8
Sk = m.Array(m.Var,nums_k)
for i in range(nums_k):
    Sk[i].LOWER = 0
    Sk[i].UPPER = 1
    Sk[i].value = Sk_initial 
    Sk[i].value = Sk[i].value.value

Ik_initial = 0.1
Ik = m.Array(m.Var,nums_k)
for i in range(nums_k):
    Ik[i].LOWER = 0
    Ik[i].UPPER = 1
    Ik[i].value = Ik_initial
    Ik[i].value = Ik[i].value.value

Vk_initial = 0.1
Vk = m.Array(m.Var,nums_k)
for i in range(nums_k):
    Vk[i].LOWER = 0
    Vk[i].UPPER = 1
    Vk[i].value = Vk_initial
    Vk[i].value = Vk[i].value.value

#set the target
J1 = m.Var()
J1_initial = 0
J1.value = J1_initial
J1.value = J1.value.value
J1.LOWER = 0

J2 = m.Var()
J2_initial = 0
J2.value = J2_initial
J2.value = J2.value.value

#set the equations
eqns = []
for i,k in enumerate(degrees):
    eqns.append(Sk[i].dt() == -k*beta*Sk[i]*sum(np.array(degrees)*pk*Ik)/k_mean-Uk[i]*Sk[i]+phi*Vk[i]+gamma*Ik[i])
    eqns.append(Ik[i].dt() == k*beta*(Sk[i]+sigma*Vk[i])*sum(np.array(degrees)*pk*Ik)/k_mean-gamma*Ik[i])
    eqns.append(Vk[i].dt() == -sigma*k*beta*Vk[i]*sum(np.array(degrees)*pk*Ik)/k_mean+Uk[i]*Sk[i]-phi*Vk[i])
eqns.append(J1.dt() == sum(Ik))
eqns.append(J2.dt() == sum(Uk))

m.Equations(eqns)

# set the time
t = np.linspace(0, T, int(T/10)+1)
m.time = t
p = np.zeros(int(T/10)+1);p[-1] = 1
final=m.Param(value=p)

# initialize with simulation
m.options.IMODE=6
m.options.NODES=2
m.options.MAX_ITER=2000
m.options.SOLVER = 3
m.options.TIME_SHIFT = 0
m.options.OTOL = 1e-3
m.options.RTOL = 1e-3


m.Minimize(J1*final)
m.solve(disp=True)

and the txt file is like this.

4 1682 0.3364
5 970 0.194
6 541 0.1082
7 370 0.074
8 282 0.0564
9 197 0.0394
10 171 0.0342
11 118 0.0236
12 101 0.0202
13 73 0.0146
14 68 0.0136
15 54 0.0108
16 40 0.008
17 37 0.0074
18 39 0.0078
19 20 0.004
20 20 0.004
21 10 0.002
22 17 0.0034
23 13 0.0026
24 9 0.0018
25 13 0.0026
26 7 0.0014
27 13 0.0026
28 11 0.0022
29 7 0.0014
30 9 0.0018
31 6 0.0012
32 8 0.0016
33 3 0.0006
34 6 0.0012
35 4 0.0008
36 3 0.0006
37 9 0.0018
38 4 0.0008
39 6 0.0012
40 5 0.001
41 4 0.0008
42 2 0.0004
44 1 0.0002
46 3 0.0006
47 1 0.0002
48 2 0.0004
49 1 0.0002
50 1 0.0002
51 4 0.0008
53 2 0.0004
56 2 0.0004
57 1 0.0002
62 1 0.0002
63 1 0.0002
64 1 0.0002
68 2 0.0004
69 2 0.0004
72 1 0.0002
76 1 0.0002
77 1 0.0002
79 1 0.0002
84 1 0.0002
86 1 0.0002
91 2 0.0004
95 1 0.0002
96 1 0.0002
97 1 0.0002
102 1 0.0002
104 2 0.0004
121 1 0.0002
122 1 0.0002
147 3 0.0006
148 2 0.0004
158 1 0.0002
187 1 0.0002

score 0 · Answer 1 · answered Dec 27 '21 at 15:20

One way to speed up the calculation is to use a non-uniform time instead of the 0-1000 with 101 time steps. The uniform time may have steps that are too big right at the beginning [0,10,20...]. Alternative time steps can improve accuracy but also enables larger steps where less resolution is needed such as [0,0.1,0.2,0.5,1,2,5,10,20...900,1000]. Time steps depend on the application so you'll need to customize it for the important features of the problem.

Uniform Time Steps

# set the time
t = np.linspace(0, T, int(T/10)+1)
m.time = t

Original problem size:

 Number of state variables:          79600
 Number of total equations: -        72400
 Number of slack variables: -            0
 ---------------------------------------
 Degrees of freedom       :           7200

Alternative Time Steps

# set the time
t = [0,0.1,0.2,0.5,1,2,5,10,20,50,100,200,300,400,500,600,700,800,900,1000]
m.time = t
p = np.zeros_like(t);p[-1] = 1
final=m.Param(value=p)

 Number of state variables:          15124
 Number of total equations: -        13756
 Number of slack variables: -            0
 ---------------------------------------
 Degrees of freedom       :           1368

A few other suggestions:

Are separate differential equations needed for each row of the data file? If the objective is to align the model with data then consider using IMODE=5. See Example 3.
There is an SEIR model with dynamic optimization (IMODE=6) here that may be helpful.
Statements such as J2.value = J2.value.value aren't needed. The values are correctly initialized with J2.value = J2_initial.
Open the run directory with m.open_folder() to view the model file gk0_model.apm.
Use the gekko summation with m.sum() instead of sum() for improved performance.
Try initialization with m.options.COLDSTART=1 to verify zero DOF for simulation and provide a good initial guess for the solver.
A reduced data set can help with application development. Try a shorter data file as the application is under development such as:

4 1682 0.3364
5 970 0.194
6 541 0.1082
7 370 0.074
8 282 0.0564
9 197 0.0394
10 171 0.0342
11 118 0.0236

Switch to a different solver such as m.options.SOLVER=1
Start with a shorter time horizon and add time points. Time steps that are too large can also cause convergence failure.

Here is a modified version that solves successfully:

import numpy as np
from gekko import GEKKO
import matplotlib.pyplot as plt
import os

N = 5000
T = 1000
beta = 0.025
phi = 0.005
gamma = 0.01
sigma = 0.001
kmin = 4
kmax = 187
k_mean = 3.9968

data = np.loadtxt("data.txt")
pk = data[:,2].tolist()
degrees = data[:,0].tolist()
nums_k = len(degrees)

m = GEKKO(remote=True)

#set the variables
Uk = m.Array(m.MV,nums_k)
for i in range(nums_k):
    Uk[i].LOWER = 0
    Uk[i].UPPER = 1
    Uk[i].STATUS = 1

Sk_initial = 0.8
Sk = m.Array(m.Var,nums_k)
for i in range(nums_k):
    Sk[i].LOWER = 0
    Sk[i].UPPER = 1
    Sk[i].value = Sk_initial 

Ik_initial = 0.1
Ik = m.Array(m.Var,nums_k)
for i in range(nums_k):
    Ik[i].LOWER = 0
    Ik[i].UPPER = 1
    Ik[i].value = Ik_initial

Vk_initial = 0.1
Vk = m.Array(m.Var,nums_k)
for i in range(nums_k):
    Vk[i].LOWER = 0
    Vk[i].UPPER = 1
    Vk[i].value = Vk_initial

#set the target
J1 = m.Var()
J1_initial = 0
J1.value = J1_initial
J1.LOWER = 0

J2 = m.Var()
J2_initial = 0
J2.value = J2_initial

#set the equations
eqns = []
for i,k in enumerate(degrees):
    eqns.append(Sk[i].dt() == -k*beta*Sk[i]*\
        m.sum(np.array(degrees)*pk*Ik)/k_mean\
            -Uk[i]*Sk[i]+phi*Vk[i]+gamma*Ik[i])
    eqns.append(Ik[i].dt() == k*beta*(Sk[i]+sigma*Vk[i])\
        *m.sum(np.array(degrees)*pk*Ik)/k_mean-gamma*Ik[i])
    eqns.append(Vk[i].dt() == -sigma*k*beta*Vk[i]\
        *m.sum(np.array(degrees)*pk*Ik)/k_mean\
            +Uk[i]*Sk[i]-phi*Vk[i])
eqns.append(J1.dt() == m.sum(Ik))
eqns.append(J2.dt() == m.sum(Uk))

m.Equations(eqns)

# set the time
#t = [0,0.1,0.2,0.5,1,2,5,10,20,50,100,200,300,400,500,600,700,800,900,1000]
t = [0,0.001,0.01,0.1,1,5,10,20,30,40,50,60,70,80,90,100]
m.time = t
p = np.zeros_like(t);p[-1] = 1
final=m.Param(value=p)

# initialize with simulation
m.options.IMODE=6
m.options.NODES=3
m.options.MAX_ITER=2000
m.options.SOLVER = 3
m.options.TIME_SHIFT = 0
m.options.OTOL = 1e-3
m.options.RTOL = 1e-3

m.Minimize(J1*final)

#m.open_folder()
m.options.SOLVER=3

# Optional: Try Coldstart first
m.options.COLDSTART=1
m.solve(disp=True)

# If Coldstart is successful then optimize
m.options.COLDSTART=0
m.solve(disp=True)

This gives a successful solution with the smaller data file:

EXIT: Optimal Solution Found.
 
 The solution was found.
 
 The final value of the objective function is    64.4192916072581     
 
 ---------------------------------------------------
 Solver         :  IPOPT (v3.12)
 Solution time  :    8.59049999999115      sec
 Objective      :    64.4192916072581     
 Successful solution
 ---------------------------------------------------

I use gekko for optimal control of complex network, but if often can't extract hessian matrix (stacked in this step)

1 Answers1