Parallelization of an optimization problem in GEKKO

user_8281509 0
Viewed 320

I'm simulating an optimization problem in GEKKO using the following code.

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

In the code, the variable tm = np.linspace(0, tf, nt) is modified to check how tm changes the estimated parameter values. When nt is more, the time taken by the solver to converge to a solution is large. So I am trying to parallelize the code. I had a look at the GEKKO example available in this tutorial. I'd like to adapt the procedure given in the above-mentioned link.

But, I could understand a few steps. For instance, in the following code provided in the link:

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

Here, all the variables are appended with s.m. . 1. Should I also append all the variables with s.m? Or just the lines that have m.something? 2. ai, bi are passed to def _init in the above code, in my example should I pass tm?

Clarifications on these doubts and explanations on how to proceed will be of great help.

EDIT: From the explanations given below and also from table 3 presented in the below-mentioned reference, I understand I should be using 1. COLDSTART=2 when the solver setting is IPOPT or 2.

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

I've tried to implement the second strategy(2). The code is not complete yet.

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

-- here, I'd like to confirm if ìnitial guess refers to guesses for parameters Dhat0 = 5000*np.ones(ngrid-1) in my code or the initial conditions of the state variables in the ode given in m.Equations.

I tried,
m.options.IMODE = 5
m.solve()
print(Dhat)

prints all 5000 which is same as input.

Further suggestions, please.

# 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)
nonlinear-optimization gekko parallel-processing python-3.x multithreading
1 Answers

There are two different types of parallel computing available in Gekko.

  • Parallel linear solvers in IPOPT with ma77, ma97, and others. This is typically only a 20-60% improvement in speed from some of the testing that I've done on large-scale problems. These options aren't available in the IPOPT version that is distributed publicly because the solvers require a license. The linear solver MUMPS is distributed with Gekko but does not include parallel support (although this is potentially coming later). The issue is that the solver is only part of the solution and even if the solver were infinitely fast, the automatic differentiation, objective evaluation, and equation evaluation still takes about 50% of the CPU time.
  • The other method for parallelization is when you have separate simulations that can run independently. This is often called "massively parallel" because the processes can be split into separate threads and then the code combines again when all of the sub-processes are complete. The link that you found uses multi-threading. Your problem is not set up for multi-threading.

If you want to improve the speed, I recommend that you try initialization with IMODE=7 and then feed that simulation solution as an initial guess for IMODE=5. An alternative is to use COLDSTART=2 and then solve the optimization problem as a next solve with COLDSTART=0 and TIME_SHIFT=0. These strategies are discussed in:

Response to Edit

Try inserting the following instead of a single m.solve() command:

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
user_23669410
Related