python - 使用 GEKKO 模拟具有巨大阵列的状态空间方程

Question

我正在尝试使用 gekko 模拟具有输入 (uu_DF) 和初始条件 (init_cond_DF) 的状态空间方程。该代码适用于小型数组(例如 3*3)。然而，在高排名的 A、B、C 矩阵(等级(A)=4553)中，代码出现以下错误请注意，我正在使用 google colab 来运行此代码，因为我的笔记本电脑(i5 CPU @ 2.67 GHz + 4GB RAM)无法运行它!

apm 35.196.244.129_gk_model5 <br><pre> ----------------------------------------------------------------
 APMonitor, Version 0.9.2
 APMonitor Optimization Suite
 ----------------------------------------------------------------
 
 
 --------- APM Model Size ------------
 Each time step contains
   Objects      :            1
   Constants    :            0
   Variables    :        27318
   Intermediates:            0
   Connections  :        27318
   Equations    :            0
   Residuals    :            0
 
 @error: Model Undefined
 Error, did not find at least one of A, B, or C matrices: 
 statespace1.a.txt
 statespace1.b.txt
 statespace1.c.txt
 Stopping...

---------------------------------------------------------------------------

Exception                                 Traceback (most recent call last)

<ipython-input-13-bb51a821614c> in <module>()
     14 #m.options.nodes = 2
     15 m.open_folder
---> 16 m.solve() # (GUI=True)
     17 
     18 

/usr/local/lib/python3.6/dist-packages/gekko/gekko.py in solve(self, disp, debug, GUI, **kwargs)
   2172             #print APM error message and die
   2173             if (debug >= 1) and ('@error' in response):
-> 2174                 raise Exception(response)
   2175 
   2176             #load results

Exception:  @error: Model Undefined
 Error, did not find at least one of A, B, or C matrices: 
 statespace1.a.txt
 statespace1.b.txt
 statespace1.c.txt
 Stopping...

我的代码如下所示

import numpy as np
from gekko import GEKKO
import matplotlib.pyplot as plt
import pandas as pd


# Importing matrices via .csv file
M = pd.read_csv(r'M.csv',header=None)
C = pd.read_csv (r'C.csv',header=None)
K = pd.read_csv (r'K.csv',header=None)
X0 = pd.read_csv (r'X0.csv',header=None)
F = pd.read_csv (r'F.csv',header=None)

# Some prerequisite calculations
dim = len(X0)
DF_eye = pd.DataFrame(np.eye(dim,dtype='uint8'))
DF_zeros = pd.DataFrame(np.zeros((dim, dim),dtype='uint8'))
BB = DF_eye # coefficient inside B matrix
M_inv_DF =  pd.DataFrame(np.linalg.inv(M))
CC = pd.DataFrame(np.zeros((2*dim, 2*dim),dtype='uint8')) # controller matrix

# System dynamics matrix
A_Row1_DF = pd.concat((DF_zeros, DF_eye), axis=1)
A_Row2_DF = pd.concat((M_inv_DF * K, M_inv_DF * C), axis=1)
A_DF = pd.concat((A_Row1_DF, A_Row2_DF), axis=0)

# System input matrix
B_Row1_DF = pd.concat((DF_zeros, DF_zeros), axis=1)
B_Row2_DF = pd.concat((DF_zeros, M_inv_DF * BB), axis=1)
B_DF = pd.concat((B_Row1_DF, B_Row2_DF), axis=0)

# Input matrix
uu_DF = pd.concat((pd.DataFrame(np.zeros((dim, 1))), pd.DataFrame(F)), axis=0)

# Initial condition matrix
init_cond_DF = pd.concat((X0, X0), axis=0)

    
m = GEKKO()  # create GEKKO model

x,y,u = m.state_space(A_DF,B_DF,CC,D=None)


m.time = np.linspace(0,120,2*dim)
u[0].value = np.array(uu_DF)
x[0].value = np.array(init_cond_DF)
m.options.imode = 4

m.solve(GUI=True) # (GUI=True)


plt.subplot(2,2,1)
plt.plot(m.time,u[0].value,'r-',label=r'$Input$')
plt.legend()
plt.subplot(2,2,2)
plt.plot(m.time,x[0].value,'b--',label=r'$Initial condition$')
plt.legend()
plt.subplot(2,2,3)
plt.plot(m.time,x[4],'.-',label=r'$x_4$')
plt.legend()
plt.show()

我针对秩的 A 矩阵测试了我的代码:82、126、170、214、258、298、342、386、430、474 前四个系统以120 秒 的最终时间 求解，求解时间 最多为4.36 秒。接下来的 五个 系统以 60 秒 的最终时间 解决，解决时间 呈指数增长到 100.13 秒。秩为 474 的 A 矩阵系统是 gekko 提供其模拟的最后一个系统，最终时间为 10 秒，解决时间 113.5 秒。 执行时间为8700.554s!在我下一次尝试902 最终时间 为 10 秒 的 A 矩阵时，gekko 引发了以下错误:

Error: 'results.json' not found. Check above for additional error details
---------------------------------------------------------------------------
FileNotFoundError                         Traceback (most recent call last)
<ipython-input-14-897bb2f533e0> in <module>()
     10 #m.options.NODES = 4
     11 #m.open_folder()
---> 12 m.solve() # (GUI=True)

________________________________________
1 frames
________________________________________
/usr/local/lib/python3.6/dist-packages/gekko/gk_post_solve.py in load_JSON(self)
     11 ## options.JSON has all APM options
     12 def load_JSON(self):
---> 13     f = open(os.path.join(self._path,'options.json'))
     14     data = json.load(f)
     15     f.close()

FileNotFoundError: [Errno 2] No such file or directory: '/tmp/tmpajn2dnt_gk_model1/options.json'

Answer 1

如果您只需要模拟状态空间模型，那么 scipy.signal 可能是一个不错的选择。 Scipy.signal here 还有其他示例.

import numpy as np
from scipy import signal

A = [[0.0,1.0],[-0.25,-0.5]]
B = [[0.0],[0.75]]
C = [[1.0,0.0]]
D = [0.0]
sys3 = signal.StateSpace(A,B,C,D)
t = np.linspace(0,30,200)
u = np.zeros(len(t))
u[5:100] = 1.0 # first step input
u[100:] = 2.0  # second step input
t3,y3,x3 = signal.lsim(sys3,u,t)

from gekko import GEKKO
m = GEKKO(remote=False)
xgk,ygk,ugk = m.state_space(A,B,C,D=None)
m.time = t; ugk[0].value = u
m.options.imode = 4
m.options.nodes = 4
m.solve(disp=False)

import matplotlib.pyplot as plt
plt.plot(t3,y3,'k-',lw=2)
plt.plot(t,ygk[0].value,'y--',lw=2)
plt.plot(t,u,'r-')
plt.legend(['y Scipy','y Gekko','u'],loc='best')
plt.xlabel('Time')
plt.show()

如果您需要使用 Gekko 进行模型预测控制，那么您可以使用 dense=False matrix option将矩阵转换为稀疏形式。这有助于大型系统。还有 Time Series (ARX)和 Discrete State Space形式。