using GEKKO to simulate a State Space equation with huge arrays

Question

I'm trying to use gekko to simulate a state space equation with input (uu_DF) and initial condition (init_cond_DF). The code works fine for small arrays (e.g., 3*3). However in high ranked A,B,C matrices (rank(A)=4553) the code comes to the following error Please note that I'm using google colab to run this code since my laptop (i5 CPU @ 2.67 GHz + 4GB RAM) is not able to run it!

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

My code is shown below

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

I tested my code for the A matrices of ranks: 82, 126, 170, 214, 258, 298, 342, 386, 430, 474 the first four systems were solved with a final time of 120 sec and the solution time was at most 4.36 sec. The next five systems were solved with a final time of 60 sec and the solution time increases exponentially to 100.13 sec. The system with the A matrix of rank 474 is the last one that gekko affords its simulation with a final time of 10 sec and the solution time of 113.5 sec. with the Execution time of 8700.554s! In my next try for the A matrix of rank 902 with a final time of 10 sec gekko raises the following error:

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

If you just need to simulate a state space model then scipy.signal may be a good option. There are additional examples with 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()

If you need to use Gekko for Model Predictive Control then you can use the dense=False matrix option that converts the matrix to a sparse form. This helps for large-scale systems. There is also the Time Series (ARX) and Discrete State Space form.