设置最优控制问题的变量设置的正确选择是什么？

计算科学 约束优化 非线性规划 数值建模 最优控制 壁虎
2021-12-21 15:33:58

这是我之前的问题的后续行动

我有以下动力系统，

(1)dϕdt=−MTDMϕ

(2)dϕ^dt=−MTD~Mϕ^

(1)表示系统的精确动力学，并且(2)是应该给出相同时间课程概况的近似动态(1), 优化后。理想情况下，我正在解决同一系统的动力学问题(1)和(2).(2)更像是一个不安的版本(1). 扰动是通过设置完成的D^= D/10。为了理解起见，让我们假设(1)给出实验值和(2)是预测值。

我已经在 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)

RESULT： Dhat是求解器返回的参数向量。Dhat由优化器拟合以最小化状态变量的测量值和预测值之间的误差。

为了检查求解器的执行情况，我设置D~（在等式 2 中，模型系统）=D（在等式 1 中，实际系统）进行初步测试。这意味着等式 1 等于等式 2（无扰动）；目标误差为零；的输出D~求解器返回的预期等于输入D，在等式 1 中。

但是，Dhat求解器返回的输出只有在Dhat代码中初始化为操纵变量 (m.MV) 而不是固定变量 (m.FV) 时才等于 D。

什么时候， Dhat = [m.MV(value=Dhat0[i]) for i in range(0, ngrid-1)]

最后一个时间点的输出：

4965.7481122
4969.8889601
4977.2097334
4991.4733925
5003.2160715   
5008.6109002
5008.2076146
5004.688907
5000.8233427
Objective      :  2.377072623938945

什么时候， Dhat = [m.FV(value=Dhat0[i]) for i in range(0, ngrid-1)]

所有时间点的输出：

3841.8094003
4198.623965
5319.3033474
6065.5329592
6467.5482342
6703.7146752
6859.9707626
9454.6282272
5098.1687634

 Objective      :  0.3068466339064452

我不确定为什么这些设置返回的解决方案存在差异以及为什么求解器不返回D~=D（在代码中为 Dhat 设置的初始值Dhat0 = 5000*np.ones(ngrid-1)）当等式 1 = 等式 2 时。

任何解释都会非常有帮助。

编辑：我也想了解搭配时间点在解决这个最优控制问题中的作用

我将时间点的数量nt从 51nt = int(tf / 0.01) + 1 更改为 501 nt = int(tf / 0.001)+ 1，但没有找到解决方案。在这里，我正在尝试检查增加是否nt会在使用时返回所有 5000m.FV

收敛失败的最后一次迭代的目标函数

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  60 6.8984929e+002 8.58e+002 2.45e+013 -11.0 2.77e+002   5.6 1.00e+000 5.00e-001h  2

WARNING: Problem in step computation; switching to emergency mode.
  63r7.3217465e+002 1.86e+002 9.99e+002   0.3 0.00e+000   6.0 0.00e+000 0.00e+000R  1

MUMPS returned INFO(1) =-13 - out of memory when trying to allocate 128655080 bytes.
In some cases it helps to decrease the value of the option "mumps_mem_percent".
WARNING: Problem in step computation; switching to emergency mode.
Restoration phase is called at point that is almost feasible,
  with constraint violation 0.000000e+000. Abort.
Restoration phase in the restoration phase failed.

从这里的报道来看，我知道mumps_mem_percent它与 IPOPT 求解器有关，但我不确定如何更改设置。我想知道如何增加mumps_mem_percentGEKKO。

编辑2：根据下面的解释，我尝试检查由集成求解器生成的解决GEKKO方案scipy's odeint

在此处输入图像描述

我可以观察到，在初始时间步长处，使用 GEKKO 中的积分求解器生成的解会产生负值。如果降低相对/绝对公差会有所帮助吗？我不确定这里使用的默认值。在此处提供的文档中，对于 scipy 的 odeint，rtol 和 atol 默认为 = 1.49012e-8。

EDIT3：按照下面的建议更改 rtol 和 otol 后，GEKKO 仍然在初始时间步返回负值。以下代码仅用于求解和比较 odeint 和 GEKKO 中的微分方程。请注意：m.options.NODES = 3 不用于仅求解和比较颂歌。

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
    """
    # M^T
    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.05
    nt = int(tf / 0.0001) + 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
    # ------------------------------------------------------------------------------------------------------------------
    m = GEKKO(remote=False)
    m.time = tm

    # ------------------------------------------------------------------------------------------------------------------
    # initialize phi_hat
    # ------------------------------------------------------------------------------------------------------------------

    phi_hat = [m.Var(value=phi_0[i]) for i in range(ngrid)]

    # ------------------------------------------------------------------------------------------------------------------
    # state variables
    # ------------------------------------------------------------------------------------------------------------------

    #phi_hat = m.CV(value=phi)
    #phi_hat.FSTATUS = 1  # fit to measurement phi obtained from 'def actual'

    # ------------------------------------------------------------------------------------------------------------------
    # parameters (/control variables to be optimized by GEKKO)
    # ref: http://apmonitor.com/do/index.php/Main/DynamicEstimation
    # def model
    # ------------------------------------------------------------------------------------------------------------------

    Dhat0 = 5000*np.ones(ngrid-1)
    Dhat = [m.FV(value=Dhat0[i]) for i in range(0, ngrid-1)]
    # Dhat.STATUS = 1  # adjustable parameter

    # ------------------------------------------------------------------------------------------------------------------
    # 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
    pprint(int_value.shape)
    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]))

    # ------------------------------------------------------------------------------------------------------------------
    # objective
    # ------------------------------------------------------------------------------------------------------------------
    # f = sum((phi(:) - phi_tilde(:)).^2);(MATLAB)
    # m.Minimize()

    # ------------------------------------------------------------------------------------------------------------------
    # simulation
    # ------------------------------------------------------------------------------------------------------------------
    m.options.IMODE = 4  # simultaneous dynamic estimation
    #m.options.NODES = 3  # collocation nodes
    #m.options.EV_TYPE = 2  # squared-error :minimize model prediction to measurement
    m.options.RTOL = 1.0e-8
    m.options.OTOL = 1.0e-8
    m.solve()

    """
    #------------------------------------------------------------------------------------------------------------------
    plt.figure()
    for i in range(0, ngrid):
        plt.plot(tm*60, phi_hat[i].value, '--', label=f'gekko_{i}')
        plt.plot(tm*60, phi[:, i], label=f'odeint_{i}')
    plt.legend(loc="upper right")
    plt.ylabel('phi/phi_hat')
    plt.xlabel('Time (s)')
    plt.xlim([0, 3])
    plt.show()

在此处输入图像描述

1个回答

该m.MV()类型具有其他调整参数，例如可能导致解决方案差异的移动抑制。此外，m.MV()在每个时间点都可以调整，而不是在整个时间窗口中调整m.time单个值。通过对 进行以下调整，m.FV()您可以获得与 类似的结果。FVMV

    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].MV_STEP_HOR = nt+10
  • 设置DCOST为零
        Dhat[i].DCOST = 0

通过这些更改，它提供了一个更接近于使用FV. 下一个问题是参数值是否应该全部5000而不是FV报告值。参数值的变化可能是由于 Gekko 和 ODEINT 求解方程的方式不同。您可以使用m.options.SENSITIVITY=1进行测试的参数敏感性可能较低（对于预测结果的微小变化，Dhat 会移动很多） 。两种解决方案也可能存在一些集成容差或离散化差异。在这些情况下，您通过移动参数来匹配两个响应以实现对齐。

对编辑 1 的回应：

本课程材料中还有关于搭配节点的更多信息。如果您在课程中搜索搭配，还会出现其他问题。为了摆脱逆响应，我建议您要么添加一个约束，>=0要么使用m.options.NODES=2并增加时间点的数量。

对编辑 2 的响应：您可以使用以下方法调整客观容差：

m.options.otol = 1e-8

以及方程解的容差：

m.options.rtol = 1e-8

在大多数情况下，从默认值调整这些容差通常没有多大帮助。如果您需要更准确的解决方案，我建议您调整问题中的节点数。增加节点的数量还具有增加解决方案的求解时间和内存需求的效果。

对编辑 3 的回应：

一开始您可能需要更多的时间点来提高准确性。有时我会创建一个带有一些 logspace 和 linspace 元素的时间向量：

import numpy as np
t1 = np.linspace(0,1,11)
print('Linear Sequence')
print(t1)
t2 = np.logspace(-3,-1.01,base=10)
print('Log10 Sequence')
print(t2)
t = np.insert(t1,1,t2)
print('Combined')
print(t)

Linear Sequence
[0.  0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1. ]
Log10 Sequence
[0.001      0.00109803 0.00120566 0.00132384 0.00145361 0.0015961
 0.00175256 0.00192436 0.00211299 0.00232012 0.00254755 0.00279727
 0.00307147 0.00337256 0.00370315 0.00406615 0.00446474 0.00490239
 0.00538295 0.00591061 0.00649    0.00712619 0.00782473 0.00859175
 0.00943396 0.01035872 0.01137413 0.01248908 0.01371333 0.01505758
 0.0165336  0.0181543  0.01993388 0.0218879  0.02403346 0.02638934
 0.02897616 0.03181655 0.03493537 0.03835991 0.04212014 0.04624897
 0.05078252 0.05576048 0.06122641 0.06722813 0.07381817 0.0810542
 0.08899954 0.09772372]
Combined
[0.         0.001      0.00109803 0.00120566 0.00132384 0.00145361
 0.0015961  0.00175256 0.00192436 0.00211299 0.00232012 0.00254755
 0.00279727 0.00307147 0.00337256 0.00370315 0.00406615 0.00446474
 0.00490239 0.00538295 0.00591061 0.00649    0.00712619 0.00782473
 0.00859175 0.00943396 0.01035872 0.01137413 0.01248908 0.01371333
 0.01505758 0.0165336  0.0181543  0.01993388 0.0218879  0.02403346
 0.02638934 0.02897616 0.03181655 0.03493537 0.03835991 0.04212014
 0.04624897 0.05078252 0.05576048 0.06122641 0.06722813 0.07381817
 0.0810542  0.08899954 0.09772372 0.1        0.2        0.3
 0.4        0.5        0.6        0.7        0.8        0.9
 1.        ]

您可能还想检查初始条件以查看是否导致负值。

其它你可能感兴趣的问题
上一篇为什么在共轭梯度中用作预处理器的平滑聚合多重网格方法会减慢求解时间？ 下一篇计算从向量到向量矩阵的距离