連立微分方程式の数値解析

制御入力のあるシステム

制御入力$u(t)$がある線形システムの挙動をシミュレーションすることを考える。
このとき、式（１）の連立微分方程式を解けばよい。

\dfrac{dx}{dt} = Ax+Bu \\

y=Cx
\tag{1}

実装

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

def func(x,t,A,B,u_list,t_list):
    x = np.array(x)
    x = x.reshape(x.shape[0],1)
    
    u = np.zeros_like(x)
    t_indx = np.where(t_list<=t)
    t_indx = np.max(t_indx[0])
    for i in range(u.shape[0]):
        u[i] = u_list[i][t_indx]
    
    A = np.array(A)
    B = np.array(B)
    
    dxdt = np.zeros_like(x)

    dxdt = A@x+B@u

    dxdt = dxdt.reshape(1,dxdt.shape[0])
    dxdt = dxdt.tolist()[0]
    return dxdt

def scenario_1():
    x_init = [10,5,2]
    t_list = np.linspace(0.0,10,1000)
    u = np.sin(np.linspace(0.0,10,1000))
    u_list  = np.array([u,u,u])

    A = [[-1,0,0],[0,-1,0],[0,0,-1]]
    B = [[1,0,0],[0,1,0],[0,0,1]]
    C = [[1,0,0],[0,1,0],[0,0,1]]

    x = odeint(func,x_init,t_list,args=(A,B,u_list,t_list))
    y = np.array(C)@x.T
    y = y.T
    return x,y,u_list

def savefig(sol,pngfn):
    plt.figure()
    for i in range(sol.shape[1]):
        plt.plot(sol[:,i])
    plt.savefig(pngfn)
    plt.close()

def savecsv(x,y,u_list,csvfn):
    sol = np.concatenate([x,y,u_list.T],axis=1)
    np.savetxt(csvfn, sol, delimiter=',')

x,y,u_list = scenario_1()
pngfn = "scenario_1.png"
savefig(y,pngfn)
csvfn = "scenario_1.csv"
savecsv(x,y,u_list,csvfn)