【原】使用Python-Gekko庫(kù)求解一階ODE

基算仿真 2023-05-30 發(fā)布于江蘇

展開(kāi)全文

之前在“使用Python求解微分方程（1）一階ODE”介紹了使用Python的Scipy庫(kù)求解一階常系數(shù)微分方程,。

今天介紹一個(gè)名為Gekko的第三方庫(kù),，以及如何使用Gekko求解一階常系數(shù)微分方程,。

Gekko介紹

Gekko(壁虎)是一個(gè)用于機(jī)器學(xué)習(xí)和動(dòng)態(tài)優(yōu)化的面向?qū)ο蟮腜ython第三方庫(kù),，它可以與線(xiàn)性,、二次,、非線(xiàn)性和混合整數(shù)編程（LP,、QP,、NLP,、MILP、MINLP）的大規(guī)模求解器相配合,?？梢詫?shí)現(xiàn)的功能包括參數(shù)回歸、數(shù)據(jù)調(diào)節(jié),、實(shí)時(shí)優(yōu)化,、動(dòng)態(tài)模擬和非線(xiàn)性預(yù)測(cè)控制。

下面介紹一些Gekko的常規(guī)操作,，方便理解后文的方程求解,。

使用Gekko時(shí)一般需要?jiǎng)?chuàng)建一個(gè)模型對(duì)象：

from gekko import GEKKOm = GEKKO([server], [name]):

定義常量

c =  m.Const(value, [name])#value必須提供且必須是標(biāo)量

定義變量

v = m.Var([value], [lb], [ub], [integer], [name])#lb 和 ub 分別為求解器提供變量下限和上限#integer 是一個(gè)布爾值，它為混合整數(shù)求解器指定一個(gè)整數(shù)變量

定義方程

m.Equation(equation)#例如m.Equation(3*x == (y**2)/z)

方程均使用隱式求解

求解案例

這里以下述微分方程為例進(jìn)行求解,。

我們令k=0.3,，y的初始值y0=5進(jìn)行計(jì)算。

求解代碼為：

import numpy as npfrom gekko import GEKKOimport matplotlib.pyplot as plt
m = GEKKO()    # create GEKKO modelk = 0.3        # constanty = m.Var(5.0) # create GEKKO variablem.Equation(y.dt()==-k*y) # create GEKKO equationm.time = np.linspace(0,20) # time points
# solve ODEm.options.IMODE = 4m.solve(disp=False)
# plot resultsplt.plot(m.time,y)plt.xlabel('time')plt.ylabel('y(t)')plt.show()

運(yùn)行結(jié)果為：

接下來(lái)可以嘗試修改系數(shù)k,，看看結(jié)果有何不同,。這里令k=0.1、0.5,、0.7分別計(jì)算試試,，求解代碼為：

import numpy as npfrom gekko import GEKKOimport matplotlib.pyplot as plt
m = GEKKO()    # create GEKKO modelk = m.Param()  # constanty = m.Var(5.0) # create GEKKO variablem.Equation(y.dt()==-k*y) # create GEKKO equationm.time = np.linspace(0,20) # time points
# solve ODEs and plotm.options.IMODE = 4m.options.TIME_SHIFT=0
k.value = 0.1m.solve(disp=False)plt.plot(m.time,y,'r-',linewidth=2,label='k=0.1')
k.value = 0.5m.solve(disp=False)plt.plot(m.time,y,'b--',linewidth=2,label='k=0.5')
k.value = 0.7m.solve(disp=False)plt.plot(m.time,y,'g:',linewidth=2,label='k=0.7')
plt.xlabel('time')plt.ylabel('y(t)')plt.legend()plt.show()

計(jì)算結(jié)果為：

再?lài)L試求解一個(gè)：

這里令u(t)為一階躍函數(shù)，在t=10時(shí)u=1,。

求解代碼為：

import numpy as npfrom gekko import GEKKOimport matplotlib.pyplot as plt
m = GEKKO()    # create GEKKO modelm.time = np.linspace(0,40,401) # time points
# create GEKKO parameter (step 0 to 1 at t=10)u_step = np.zeros(401)u_step[100:] = 2.0u = m.Param(value=u_step)  
y = m.Var(1.0) # create GEKKO variablem.Equation(5 * y.dt()==-y+u) # create GEKKO equation
# solve ODEm.options.IMODE = 4m.solve(disp=False)
# plot resultsplt.plot(m.time,y,'r-',label='Output (y(t))')plt.plot(m.time,u,'b-',label='Input (u(t))')plt.ylabel('values')plt.xlabel('time')plt.legend(loc='best')plt.show()