GEKKO中一个优化问题的并行化

Question

我正在使用以下代码模拟GEKKO中的优化问题。

# Copyright 2020, Natasha, All rights reserved.
import numpy as np
from gekko import GEKKO
from pprint import pprint
import matplotlib.pyplot as plt
from scipy.integrate import odeint

def get_mmt():
"""
M and M transpose required for differential equations
:params: None
:return: M transpose and M -- 2D arrays ~ matrices
"""
MT = np.array([[-1, 0, 0, 0, 0, 0, 0, 0, 0],
[1, -1, 0, 0, 0, 0, 0, 0, 0],
[0, 1, -1, 0, 0, 0, 0, 0, 0],
[0, 0, 1, -1, 0, 0, 0, 0, 0],
[0, 0, 0, 1, -1, 0, 0, 0, 0],
[0, 0, 0, 0, 1, -1, 0, 0, 0],
[0, 0, 0, 0, 0, 1, -1, 0, 0],
[0, 0, 0, 0, 0, 0, 1, -1, 0],
[0, 0, 0, 0, 0, 0, 0, 1, -1],
[0, 0, 0, 0, 0, 0, 0, 0, 1]])
M = np.transpose(MT)
return M, MT

def actual(phi, t):
"""
Actual system/ Experimental measures
:param  phi: 1D array
:return: time course of variable phi -- 2D arrays ~ matrices
"""
# spatial nodes
ngrid = 10
end = -1
M, MT = get_mmt()
D = 5000*np.ones(ngrid-1)
A = [e-mail](D)@M
A = A[1:ngrid-1]
# differential equations
dphi = np.zeros(ngrid)
# first node
dphi[0] = 0
# interior nodes
dphi[1:end] = -A@phi  # value at interior nodes
# terminal node
dphi[end] = D[end]*2*(phi[end-1] - phi[end])
return dphi

if __name__ == '__main__':
# ref: https://apmonitor.com/do/index.php/Main/PartialDifferentialEquations
ngrid = 10  # spatial discretization
end = -1
# integrator settings (for ode solver)
tf = 0.5
nt = int(tf / 0.01) + 1
tm = np.linspace(0, tf, nt)
# ------------------------------------------------------------------------------------------------------------------
# measurements
# ref: https://www.youtube.com/watch?v=xOzjeBaNfgo
# using odeint to solve the differential equations of the actual system
# ------------------------------------------------------------------------------------------------------------------
phi_0 = np.array([5, 0, 0, 0, 0, 0, 0, 0, 0, 0])
phi = odeint(actual, phi_0, tm)
# plot results
plt.figure()
plt.plot(tm*60, phi[:, :])
plt.ylabel('phi')
plt.xlabel('Time (s)')
plt.show()
# ------------------------------------------------------------------------------------------------------------------
#  GEKKO model
# ------------------------------------------------------------------------------------------------------------------
m = GEKKO(remote=False)
m.time = tm
# ------------------------------------------------------------------------------------------------------------------
# initialize state variables: phi_hat
# ref: https://apmonitor.com/do/uploads/Main/estimate_hiv.zip
# ------------------------------------------------------------------------------------------------------------------
phi_hat = [m.CV(value=phi_0[i]) for i in range(ngrid)]  # initialize phi_hat; variable to match with measurement
# ------------------------------------------------------------------------------------------------------------------
# parameters (/control parameters to be optimized while minimizing the cost function in GEKKO)
# ref: http://apmonitor.com/do/index.php/Main/DynamicEstimation
# ref: https://apmonitor.com/do/index.php/Main/EstimatorObjective
# def model
# ------------------------------------------------------------------------------------------------------------------
#  Manually enter guesses for parameters
Dhat0 = 5000*np.ones(ngrid-1)
Dhat = [m.MV(value=Dhat0[i]) for i in range(0, ngrid-1)]
for i in range(ngrid-1):
Dhat[i].STATUS = 1  # Allow optimizer to fit these values
# Dhat[i].LOWER = 0
# ------------------------------------------------------------------------------------------------------------------
# differential equations
# ------------------------------------------------------------------------------------------------------------------
M, MT = get_mmt()
A = MT @ np.diag(Dhat) @ M
A = A[1:ngrid - 1]
# first node
m.Equation(phi_hat[0].dt() == 0)
# interior nodes
int_value = -A @ phi_hat  # function value at interior nodes
m.Equations(phi_hat[i].dt() == int_value[i] for i in range(0, ngrid-2))
# terminal node
m.Equation(phi_hat[ngrid-1].dt() == Dhat[end] * 2 * (phi_hat[end-1] - phi_hat[end]))
# ------------------------------------------------------------------------------------------------------------------
# simulation
# ------------------------------------------------------------------------------------------------------------------
m.options.IMODE = 5  # simultaneous dynamic estimation
m.options.NODES = 3  # collocation nodes
m.options.EV_TYPE = 2  # squared-error :minimize model prediction to measurement
for i in range(ngrid):
phi_hat[i].FSTATUS = 1  # fit to measurement phi obtained from 'def actual'
phi_hat[i].STATUS = 1  # build objective function to match measurement and prediction
phi_hat[i].value = phi[:, i]
m.solve()
pprint(Dhat)

在代码中，变量tm = np.linspace(0, tf, nt)被修改以检查tm如何改变估计的参数值。当nt较大时，解算器收敛到解所花费的时间较大。所以我正在尝试将代码并行化。我看过本教程中提供的GEKKO示例。我想调整上述链接中给出的程序。

但是，我可以理解一些步骤。例如，在链接中提供的以下代码中：

def __init__(self, id, server, ai, bi):
s = self
s.id = id
s.server = server
s.m = GEKKO()
s.a = ai
s.b = bi
s.objective = float('NaN')
# initialize variables
s.m.x1 = s.m.Var(1,lb=1,ub=5)
s.m.x2 = s.m.Var(5,lb=1,ub=5)
s.m.x3 = s.m.Var(5,lb=1,ub=5)
s.m.x4 = s.m.Var(1,lb=1,ub=5)
# Equations
s.m.Equation(s.m.x1*s.m.x2*s.m.x3*s.m.x4>=s.a)
s.m.Equation(s.m.x1**2+s.m.x2**2+s.m.x3**2+s.m.x4**2==s.b)
# Objective
s.m.Obj(s.m.x1*s.m.x4*(s.m.x1+s.m.x2+s.m.x3)+s.m.x3)
# Set global options
s.m.options.IMODE = 3 # steady state optimization
s.m.options.SOLVER = 1 # APOPT solver

在这里，所有变量都附加了s.m。1.我是否也应该在所有变量后面加上s.m？或者只是有m.something的线条？2.ai、bi在上面的代码中被传递给def_init，在我的例子中，我应该传递tm吗？

澄清这些疑虑并解释如何进行将有很大帮助。

编辑：根据下面给出的解释以及下面提到的参考资料中的表3，我理解我应该使用1.当解算器设置为IPOPT或2.

IMODE=7初始化，然后提供模拟解决方案作为IMODE＝5的初始猜测。

我已经尝试实施第二种策略(2(。代码尚未完成。

将该模拟解决方案作为IMODE=5的初始猜测

这里，我想确认ìnitial guess是否指的是对参数的猜测我的代码中的Dhat0 = 5000*np.ones(ngrid-1)或m中给出的ode中状态变量的初始条件

I tried,
m.options.IMODE = 5
m.solve()
print(Dhat) 

打印与输入相同的全部5000。

请提供进一步的建议。

# Copyright 2013, Natasha, All rights reserved.
import numpy as np
from gekko import GEKKO
from pprint import pprint
import matplotlib.pyplot as plt
from scipy.integrate import odeint

def get_mmt():
"""
M and M transpose required for differential equations
:params: None
:return: M transpose and M -- 2D arrays ~ matrices
"""
MT = np.array([[-1, 0, 0, 0, 0, 0, 0, 0, 0],
[1, -1, 0, 0, 0, 0, 0, 0, 0],
[0, 1, -1, 0, 0, 0, 0, 0, 0],
[0, 0, 1, -1, 0, 0, 0, 0, 0],
[0, 0, 0, 1, -1, 0, 0, 0, 0],
[0, 0, 0, 0, 1, -1, 0, 0, 0],
[0, 0, 0, 0, 0, 1, -1, 0, 0],
[0, 0, 0, 0, 0, 0, 1, -1, 0],
[0, 0, 0, 0, 0, 0, 0, 1, -1],
[0, 0, 0, 0, 0, 0, 0, 0, 1]])
M = np.transpose(MT)
return M, MT

def actual(phi, t):
"""
Actual system/ Experimental measures
:param  phi: 1D array
:return: time course of variable phi -- 2D arrays ~ matrices
"""
# spatial nodes
ngrid = 10
end = -1
M, MT = get_mmt()
D = 5000*np.ones(ngrid-1)
A = [e-mail](D)@M
A = A[1:ngrid-1]
# differential equations
dphi = np.zeros(ngrid)
# first node
dphi[0] = 0
# interior nodes
dphi[1:end] = -A@phi  # value at interior nodes
# terminal node
dphi[end] = D[end]*2*(phi[end-1] - phi[end])
return dphi

if __name__ == '__main__':
# ref: https://apmonitor.com/do/index.php/Main/PartialDifferentialEquations
ngrid = 10  # spatial discretization
end = -1
# integrator settings (for ode solver)
tf = 0.5
nt = int(tf / 0.01) + 1
tm = np.linspace(0, tf, nt)
# ------------------------------------------------------------------------------------------------------------------
# measurements
# ref: https://www.youtube.com/watch?v=xOzjeBaNfgo
# using odeint to solve the differential equations of the actual system
# ------------------------------------------------------------------------------------------------------------------
phi_0 = np.array([5, 0, 0, 0, 0, 0, 0, 0, 0, 0])
phi = odeint(actual, phi_0, tm)
# ------------------------------------------------------------------------------------------------------------------
#  GEKKO model
# https://apmonitor.com/wiki/index.php/Main/Simulation
# ------------------------------------------------------------------------------------------------------------------
# Initialize GEKKO
m1 = GEKKO(remote=False)
m2 = GEKKO(remote=False)
for m in [m1,m2]:
m.time = tm
# ------------------------------------------------------------------------------------------------------------------
# initialize state variables: phi_hat
# ref: https://apmonitor.com/do/uploads/Main/estimate_hiv.zip
# ------------------------------------------------------------------------------------------------------------------
phi_hat = [m.CV(value=phi_0[i]) for i in range(ngrid)]  # initialize phi_hat; variable to match with measurement
# ------------------------------------------------------------------------------------------------------------------
# parameters (/control parameters to be optimized while minimizing the cost function in GEKKO)
# ref: http://apmonitor.com/do/index.php/Main/DynamicEstimation
# ref: https://apmonitor.com/do/index.php/Main/EstimatorObjective
# def model
# ------------------------------------------------------------------------------------------------------------------
#  Manually enter guesses for parameters
Dhat0 = 5000*np.ones(ngrid-1)
Dhat = [m.FV(value=Dhat0[i]) for i in range(0, ngrid-1)]
for i in range(ngrid-1):
Dhat[i].STATUS = 1  # Allow optimizer to fit these values
# Dhat[i].LOWER = 0
# ------------------------------------------------------------------------------------------------------------------
# differential equations
# ------------------------------------------------------------------------------------------------------------------
M, MT = get_mmt()
A = MT @ np.diag(Dhat) @ M
A = A[1:ngrid - 1]
# first node
m.Equation(phi_hat[0].dt() == 0)
# interior nodes
int_value = -A @ phi_hat  # function value at interior nodes
m.Equations(phi_hat[i].dt() == int_value[i] for i in range(0, ngrid-2))
# terminal node
m.Equation(phi_hat[ngrid-1].dt() == Dhat[end]*2*(phi_hat[end-1] - phi_hat[end]))
# ------------------------------------------------------------------------------------------------------------------
# simulation
# ------------------------------------------------------------------------------------------------------------------
m.options.NODES = 3  # collocation nodes
m.options.EV_TYPE = 2  # squared-error :minimize model prediction to measurement
m.options.SOLVER = 3  # 1 APOPT, 2 BPOPT, 3 IPOPT
# ------------------------------------------------------------------------------------------------------------------
#  Initialization
#  Ref: Initialization strategies for optimization of dynamic systems
# ------------------------------------------------------------------------------------------------------------------
m1.options.IMODE = 7  # simultaneous dynamic estimation
for i in range(ngrid):
phi_hat[i].FSTATUS = 1  # fit to measurement phi obtained from 'def actual'
phi_hat[i].STATUS = 1  # build objective function to match measurement and prediction
phi_hat[i].value = phi[:, i]
m1.solve()
# ------------------------------------------------------------------------------------------------------------------
#  MPH
#  Ref: Initialization strategies for optimization of dynamic systems
# ------------------------------------------------------------------------------------------------------------------
m2.options.IMODE = 5  # simultaneous dynamic estimation
for i in range(ngrid):
phi_hat[i].FSTATUS = 1  # fit to measurement phi obtained from 'def actual'
phi_hat[i].STATUS = 1  # build objective function to match measurement and prediction
phi_hat[i].value = phi[:, i]
m2.solve()
pprint(Dhat)

Answer 1

Gekko中有两种不同类型的并行计算。

• IPOPT中具有ma77、ma97等的并行线性解算器。这通常只比我在大规模问题上所做的一些测试提高了20-60%的速度。这些选项在公开发布的IPOPT版本中不可用，因为解算器需要许可证。线性求解器MUMPS与Gekko一起分布，但不包括并行支持(尽管这可能会在稍后出现(。问题是，求解器只是解决方案的一部分，即使求解器速度无限快，自动微分、客观评估和方程评估仍然需要大约50%的CPU时间
• 并行化的另一种方法是当您有可以独立运行的单独模拟时。这通常被称为"大规模并行"，因为进程可以被拆分为单独的线程，然后当所有子进程完成时，代码再次组合。您找到的链接使用多线程。您的问题没有设置为多线程

如果您想提高速度，我建议您尝试使用IMODE=7进行初始化，然后将该模拟解决方案作为IMODE=5的初始猜测。另一种选择是使用COLDSTART=2，然后将优化问题作为COLDSTART=0和TIME_SHIFT=0的下一个解决方案来解决。这些策略在以下内容中进行了讨论：

• Safdarnejad，S.M.，Hedengren，J.D.，Lewis，N.R.，Haseltine，E.，动态系统优化初始化策略，计算机和化学工程，2015，第78卷，第39-50页，DOI:10.1016/J.com.pchemeng.2015.04.016

对编辑的响应

尝试插入以下命令而不是单个m.solve()命令：

m.options.IMODE = 5      # simultaneous estimation
m.options.COLDSTART = 2  # coldstart on
m.solve(disp=False)      # first solve
m.options.COLDSTART = 0  # coldstart off
m.options.TIME_SHIFT = 0 # turn off time-shift (default=1)
m.solve(disp=False)      # second solve