python-3.x - GEKKO中优化问题的并行化

标签 python-3.x multithreading parallel-processing nonlinear-optimization gekko

我正在使用以下代码在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 = MT@np.diag(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更大时，求解器收敛到一个解所花费的时间就很大。所以我正在尝试并行化代码。我看了this教程中可用的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吗？还是只有一行东西？
2.在上面的代码中将ai，bi传递给def _init，在我的示例中，我应该传递tm吗？

澄清这些疑惑，并就如何进行进行解释，将大有帮助。

编辑:从下面给出的解释以及从下面提到的引用文献中提供的表3中，我了解我应该使用
1.当求解器设置为IPOPT或COLDSTART = 2时
2。

initialization with IMODE=7 and then feed that simulation solution as an initial guess for IMODE=5.



我已经尝试实现第二个策略(2)。代码尚未完成。

feed that simulation solution as an initial guess for IMODE=5



-在这里，我想确认ìnitial guess是否引用了对参数的猜测
我的代码中的Dhat0 = 5000*np.ones(ngrid-1)或m.Equations中给出的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 = MT@np.diag(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)

最佳答案

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.，Initialization Strategies for Optimization of Dynamic Systems，Computers and Chemical Engineering，2015，Vol。 78，第39-50页，DOI:10.1016/j.compchemeng.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

    关于python-3.x - GEKKO中优化问题的并行化，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/61143557/

    上一篇：c# - 锁定列表会产生死锁

    下一篇：javascript - 匹配的用户流Nodejs

    相关文章：

    python - Python 中迭代器的静态行为

    Java(FX) 在播放一种声音的同时播放另一种声音

    Python:为什么从线程调用的 `sys.exit(msg)` 不将 `msg` 打印到 stderr？

    Java:从同步块(synchronized block)启动新线程时会发生什么？

    cuda - 使用 CUDA 计算矩阵对应行之间的欧几里德距离

    go - 如何在 Golang 中实现适当的并行性？ goroutines 在 1.5+ 版本中是并行的吗？

    python - 为什么我的函数落后一键运行？

    python - 使用 PySpark 将日期和时间字符串转换为时间戳时如何保留毫秒？

    python - Python 倒排索引

    c - MPI 发送接收

    ©2024 IT工具网  联系我们