3
$\begingroup$

I have the following dynamical system,

$\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}$ represents the exact dynamics of a system and $\eqref{2}$ is the approximate dynamics that should give the same time course profiles as $\eqref{1}$, after optimization. Ideally, I am solving for the dynamics of the same system in $\eqref{1}$ and $\eqref{2}$. $\eqref{2}$ is more like a perturbed version of $\eqref{1}$. The perturbation is done by setting $\hat{D}$ = D/10. And for the sake of understanding, let us assume $\eqref{1}$ gives experimental values and $\eqref{2}$ are the predicted values.

The objective function includes a cost function that minimizes the difference between state variables $\phi$ and $\hat{\phi}$, by optimizing parameter $\tilde{D}$ which are the control variables.

I'm trying to solve this as a parameter estimation problem with non-linear equality constraints/defects obtained by discretizing $\eqref{2}$ at collocation points.

In MATLAB my objective function looks like the following 

[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

I've tried to set up the same problem in GEKKO. But I am not sure how to set up the objective function. [t, phi] = ode15s(@(t,phi) actual(t,phi), tspan , phi0 ,options); in MATLAB computes the time-course profiles of phi. In python code, the differential equations in function def actual(): is solved using odeint from scipy in line 102. Similarly, [t, phi_tilde] = ode15s(@(t,phi_tilde) model(t,phi_tilde, Dhat), tspan , phi0 ,options); computes the time-course profiles of phi_hat. In GEKKO, the equations of model has been set up in function def model():.

I'm stuck at this point. It's not clear to me how model that has the control variables in 1d array Dhat has to be set up and solved to compute the squared-error in loss function defined in objective function f = 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 = [email protected](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)

In short, I'd like to ask for inputs on how to set up the m.Obj, m.CV, and m.FV in GEKKO to solve this problem.

EDIT: m.CV, and m.FV have been updated in the code. I'd like to request for help with setting up the objective function alone.

Improve this question
$\endgroup$
1

1 Answer 1

Reset to default
3
$\begingroup$

A good first step with any parameter estimation problem is to solve it in simulation to verify that you can get a good solution and that the parameters have an effect on the objective. You can first simulate with m.options.IMODE=7.

Once you have an initial solution, you can set your objective function with:

for i in range(n):
    m.Minimize((phi[i]-phi_hat[i])**2)

You can have as many Minimize or Maximize statements as you need. Gekko adds all of them together to create a scalar objective value.

Alternatively, you can use the built-in objective function that comes with a CV object with m.options.EV_TYPE=2 to set a squared objective. You just need to give it values such as:

phi = m.Array(m.CV,n)
for i in range(n):
    phi[i].value = phi_hat[i]
    phi[i].FSTATUS = 1

The length of phi_hat[i] must be consistent with the number of time steps that you have defined with m.time. I recommend tutorials at https://apmonitor.com/do for more information on setting up your problem.

Improve this answer
$\endgroup$
2
  • 1
    $\begingroup$ Many thanks for the response. phi = m.Array(m.CV,n) = array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0], dtype=object) when n=10. Is it correct to understand that the value of phi is initilaized to zeros in this step. Could you please explain how to initialize it to a different value other than zeros? Doing , phi = [m.CV(value=phi_0[i]) for i in range(n)] initializes phi as a list. $\endgroup$
    – Natasha
    Commented Apr 7, 2020 at 17:10
  • 1
    $\begingroup$ You could do the following: [phi[i].value=phi_0[i] for i in range(n)]. This would only change the value and not redefine phi. $\endgroup$
    – John Hedengren
    Commented Apr 8, 2020 at 16:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.