我有以下动力系统，

$\frac{d \phi}{dt} = -M^TDM\phi \tag{1}\label{1}$

$\frac{d \hat\phi}{dt} = -M^T\tilde{D}M\hat \phi \tag{2} \label{2}$

$\eqref{1}$表示系统的精确动态，相同的时间课程配置文件的近似动态。理想情况下，我正在解决和中同一系统的动力学问题。更像是的扰动版本。扰动是通过设置 = D/10 来完成的。为了便于理解，我们假设给出了实验值，而是预测值。$\eqref{2}$$\eqref{1}$$\eqref{1}$$\eqref{2}$$\eqref{2}$$\eqref{1}$$\hat{D}$$\eqref{1}$$\eqref{2}$

目标函数包括一个成本函数，它 通过优化作为控制变量的和$\phi$$\hat{\phi}$$\tilde{D}$

我试图将此作为参数估计问题解决，该问题具有通过在搭配点$\eqref{2}$

在 MATLAB 中，我的目标函数如下所示

[Dhat,~,~,output] = fmincon(@objfun,Dhat0,[],[],[],[],[],[],@defects, opts_fmin)

function f = objfun(Dhat)

%% Integrator settings
phi0    = [5; 0; 0; 0; 0; 0; 0; 0; 0; 0];
tspan   = 0:dt:0.5;
options = odeset('abstol', 1e-10, 'reltol', 1e-9);

%% generate exact solution
    [t, phi]  = ode15s(@(t,phi) actual(t,phi), tspan , phi0 ,options);


%% generate approximate solution

    [t, phi_tilde]  = ode15s(@(t,phi_tilde) model(t,phi_tilde, Dhat), tspan , phi0 ,options);


%% objective function for fminunc/fmincon
      f = sum((phi(:) - phi_tilde(:)).^2);
end

我试图在 GEKKO 中设置同样的问题。但我不确定如何设置目标函数。[t, phi] = ode15s(@(t,phi) actual(t,phi), tspan , phi0 ,options);在 MATLAB 中计算phi. 在 python 代码中，函数def actual():中的微分方程在第 102 行使用 scipy 中的 odeint 求解。类似地，[t, phi_tilde] = ode15s(@(t,phi_tilde) model(t,phi_tilde, Dhat), tspan , phi0 ,options);计算phi_hat. 在 GEKKO 中，方程model已在 function 中建立def model():。

我被困在这一点上。我不清楚如何设置和解决一维model数组中的控制变量 以计算目标函数中定义的损失函数的平方误差。 Dhatf = sum((phi(:) - phi_tilde(:)).^2);(MATLAB)

# 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)

简而言之，我想请教如何在 GEKKO 中设置m.Obj、m.CV和m.FV来解决这个问题。

编辑： m.CV，并m.FV 已在代码中更新。我想请求帮助单独设置目标函数。