4
$\begingroup$

This is a followup to my previous question here

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.

I've set up this system in 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 = [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)

RESULT: Dhat is the parameter vector that is returned by the solver. Dhat is fit by the optimizer to minimize the error between measured and predicted values of state variables.

To check how the solver performs, I set $\tilde{D}$ (in equation 2, model system) = $D$ (in equation 1, actual system) for preliminary tests. This would imply, equation 1 is equal to equation 2 (no perturbation); error in objective will be zero; the output of $\tilde{D}$ returned by the solver will be expected to be equal to input $D$, in equation 1.

However, the output Dhat returned by the solver is equal to D only when Dhat is initialized as a manipulated variable (m.MV) in the code and not as a fixed variable (m.FV).

When, Dhat = [m.MV(value=Dhat0[i]) for i in range(0, ngrid-1)]

Output at last time point:

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

When, Dhat = [m.FV(value=Dhat0[i]) for i in range(0, ngrid-1)]

Output at all time points:

3841.8094003
4198.623965
5319.3033474
6065.5329592
6467.5482342
6703.7146752
6859.9707626
9454.6282272
5098.1687634

 Objective      :  0.3068466339064452

I am not sure why there is a difference in solution returned for these settings and why the solver doesn't return $\tilde{D} = D$ (initial values set for Dhat in the code Dhat0 = 5000*np.ones(ngrid-1)) when equation1 = equation 2.

Any explanations will be really helpful.

EDIT : I'd also like to understand the role of collocation time points in solving this optimal control problem

I changed the number of time points, nt from 51 nt = int(tf / 0.01) + 1 to 501 nt = int(tf / 0.001) + 1 and the solution was not found. Here, I am trying to check whether increasing nt will return all 5000 while using m.FV

Objective function at the last iteration where convergence failed

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.

From what's reported here, I understand mumps_mem_percent is associated with IPOPT solver, but I am not sure how to change in the settings. I'd like to know how to increase mumps_mem_percent in GEKKO.

EDIT 2: From what has been explained below, I tried to check the in solutions generated by integration solvers used in GEKKO and scipy's odeint

enter image description here

I could observe that at initial time steps solution generated by using the integration solver in GEKKO yields negative values. Would it help if the relative/ absolute tolerance is decreased? I'm not sure of the default values used here. In the documentation available here rtol and atol is = 1.49012e-8 by default for scipy's odeint.

EDIT3: After changing rtol and otol as suggested below, GEKKO still returns negative values at the initial time steps. The following code is used to solve and compare just the differential equations in odeint and GEKKO. Please note: m.options.NODES = 3 is not used to solve and compare just the odes.

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 = [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.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()

enter image description here

Improve this question
$\endgroup$

1 Answer 1

Reset to default
3
$\begingroup$

The m.MV() type has additional tuning parameters such as move suppression that is likely contributing to the difference in solution. Also, the m.MV() is adjustable at every time point in m.time instead of just a single value with an m.FV() over the entire time window. You can get similar results to an FV by making the following adjustments to the MV.

  • Set the MV_STEP_HOR to a value greater than the horizon window:
    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
  • Set the DCOST to zero
        Dhat[i].DCOST = 0

With these changes, it give a solution closer to the solution using FV. The next question is whether the parameter values should all be 5000 instead of the FV reported values. The change in parameter values is likely due to the differences in how Gekko and ODEINT solve the equations. There may be a low parameter sensitivity (Dhat moves a lot for small changes in the predicted outcome) that you can test with m.options.SENSITIVITY=1. There may also be some integration tolerance or discretization differences the two solutions. In these cases, you are matching the two responses by moving the parameters to achieve alignment.

Response to Edit 1:

There is additional information on collocation nodes with this course material. There are also additional problems if you search for collocation in the course. To get rid of the inverse response, I recommend that you either add a constraint for >=0 or else use m.options.NODES=2 and increase the number of time points.

Response to Edit 2: You can adjust the objective tolerance with:

m.options.otol = 1e-8

and the tolerance for the equation solutions with:

m.options.rtol = 1e-8

Adjusting these tolerances from defaults doesn't typically help a lot in most cases. If you need a more accurate solutions, I recommend adjusting the number of nodes that you have in your problem. Increasing the number of nodes also has the effect of increasing the solution time and memory requirements for the solution.

Response to Edit 3:

You may need a few more time points at the beginning to improve the accuracy. Sometime I create a time vector with some logspace and linspace elements:

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.        ]

You may also want to check the initial conditions to see if that is causing the negative values.

Improve this answer
$\endgroup$
7
  • $\begingroup$ Many thanks for the response. Could you please explain how the integration tolerances can be viewed in GEKKO? In the documentation available here rtol and atol is = 1.49012e-8 y default. for scipy's odeint $\endgroup$
    – Natasha
    Apr 10, 2020 at 5:40
  • $\begingroup$ I'd also like to know what additional settings can be included for m.FV to obtain results given by m.MV posted in my original post $\endgroup$
    – Natasha
    Apr 10, 2020 at 10:27
  • 1
    $\begingroup$ Many thanks for the response. I tried to simulate just the differential system, without any constraints. Please check edit. I still see negative values at initial time steps. $\endgroup$
    – Natasha
    Apr 10, 2020 at 16:44
  • 1
    $\begingroup$ Thanks a lot. You may need a few more time points at the beginning to improve the accuracy. I could see a difference in the solution curves. But could you please explain why this improves the accuracy? I've an understanding that these time points are merely used to output the values and doesn't have any effect on the time step taken by the solver in scipy's odeint and MATLAB's ode45 or 15s etc.. $\endgroup$
    – Natasha
    Apr 12, 2020 at 6:23
  • 1
    $\begingroup$ Gekko only uses dynamic time-point selection with IMODE=7 to avoid solution convergence problems. For all other cases, it only solves at the requested time points (unlike ODE solvers). This static method is a weakness and a benefit. It helps with stiff systems to automatically implement a quasi-steady-state assumption and speed up simulation or optimization convergence. It hurts with maintaining accuracy, especially at the beginning of the simulation. $\endgroup$
    – John Hedengren
    Apr 12, 2020 at 15:27

Your Answer

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

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