Распараллеливание задачи оптимизации в GEKKO - PullRequest
Новая
       64

Распараллеливание задачи оптимизации в GEKKO

2 голосов
/ 10 апреля 2020

Я моделирую задачу оптимизации в 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)

В этом коде переменная tm = np.linspace(0, tf, nt) модифицируется, чтобы проверить, как tm изменяет расчетные значения параметров. Когда nt больше, время, необходимое солверу, чтобы сходиться к решению, велико. Поэтому я пытаюсь распараллелить код. Я посмотрел пример GEKKO, доступный в этом руководстве. Я хотел бы адаптировать процедуру, приведенную в вышеупомянутой ссылке.

Но я мог бы понять несколько шагов. Например, в следующем коде, представленном в ссылке: 

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

Здесь все переменные добавляются с помощью sm. 1. Должен ли я также добавить все переменные с sm? Или только строки, которые имеют m.something? 2. ai, bi передаются в def _init в приведенном выше коде, в моем примере я должен передать tm?

Разъяснения по поводу этих сомнений и объяснения того, как действовать, будут очень полезны.

РЕДАКТИРОВАТЬ: Исходя из пояснений, приведенных ниже, а также из таблицы 3, представленной в нижеприведенной ссылке, я понимаю, что следует использовать 1. COLDSTART = 2, когда для параметра решателя установлено значение IPOPT или 2.

инициализация с IMODE = 7, а затем подать это решение для моделирования в качестве первоначального предположения для IMODE = 5.

Я пытался реализовать вторую стратегию (2). Код еще не завершен.

укажите это решение для моделирования в качестве первоначального предположения для IMODE = 5

- здесь я хотел бы подтвердить, если ìnitial guess относится к предположениям о параметрах Dhat0 = 5000*np.ones(ngrid-1) в моем коде или начальных условиях переменных состояния в ode, указанных в m.Equations. 

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

печатает все 5000, что совпадает с вводом.

Дополнительные предложения, пожалуйста. 

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

1 Ответ

2 голосов
/ 11 апреля 2020

В Gekko доступны два различных типа параллельных вычислений.

  • Параллельные линейные решатели в IPOPT с ma77, ma97 и другими. Как правило, это только на 20-60% улучшение скорости по сравнению с некоторыми тестами, которые я проводил для крупномасштабных задач. Эти параметры недоступны в версии IPOPT, распространяемой публично, поскольку решатели требуют лицензии. Линейный решатель MUMPS распространяется с Gekko, но не включает параллельную поддержку (хотя это может произойти позже). Проблема заключается в том, что решатель является лишь частью решения, и даже если решатель работает бесконечно быстро, автоматическое дифференцирование c, объективная оценка и оценка уравнений по-прежнему занимают около 50% процессорного времени.
  • Другой метод распараллеливания - это когда у вас есть отдельные симуляции, которые могут выполняться независимо. Это часто называют «массово параллельными», потому что процессы могут быть разделены на отдельные потоки, а затем код снова объединяется, когда все подпроцессы завершены. Ссылка, которую вы нашли, использует многопоточность. Ваша проблема не настроена на многопоточность.

Если вы хотите улучшить скорость, я рекомендую вам попробовать инициализацию с помощью IMODE=7, а затем передать это решение для моделирования в качестве первоначального предположения для IMODE=5. Альтернативой является использование COLDSTART=2, а затем решение задачи оптимизации в качестве следующего решения с помощью COLDSTART=0 и TIME_SHIFT=0. Эти стратегии обсуждаются в:

Ответ на редактирование

Попробуйте вставить следующее вместо одна m.solve() команда: 

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