Python：整数項を持つ常微分方程式を解く方法

Question

私はこの質問が以前に尋ねられたかどうかは分かりませんが、私は先に進んでここに投稿します。私は単純なシステムをPID私の微分方程式系は以下の通りです。私は基本的に非常に基本的なPIDアルゴリズムをコーディングしようとしています。私の制御uの構造は、誤差項の導関数と積分の両方に依存する。私は、派生語には何の問題もありません、それは私のコードで問題を作り出している積分用語です。最初に にs = 0を割り当てて、以下のコードで説明したように私の関数で使用すると、この問題は解決します。それをバイパスする方法はありますか？私はsを割り当てることを試み、グローバル変数として伝えましたが、それは私の問題を解決しませんでした。一言で言えば、私がやっているのは、毎回状態x1を追加し、dt（これはt-toldで示される）を掛けることです。

親切に私は私のコードは、以下に添付、PFAをこの問題を鉄助けます。

import numpy as np 
from scipy.integrate import ode 
import matplotlib.pyplot as plt 
plt.style.use('bmh') 

t0=0 
y0=[0.1,0.2] 
kp,kd,ki=2,0.5,0.8 
s,told=0,0 
def pid(t,Y): 
    x1,x2=Y[0],Y[1] 
    e=x1-1 
    de=x2 
    s=(x1+s) 
    integral=s*(t-told) 
    told=t 
    #ie= 
    u=kp*e+kd*de+ki*integral 
    x1dot=x2 
    x2dot=u-5*x1-2*x2 
    return[x1dot,x2dot] 

solver=ode(pid).set_integrator('dopri5',rtol=1e-6,method='bdf',nsteps=1e5,max_step=1e-3) 
solver.set_initial_value(y0,t0) 
t1=10 
dt=5e-3 
sol = [ [yy] for yy in y0 ] 
t=[t0] 
while solver.successful() and solver.t<t1: 
    solver.integrate(solver.t+dt) 
    for k in range(2): sol[k].append(solver.y[k]); 
    t.append(solver.t) 
    print(len(sol[0])) 
    print(len(t)) 
x1=np.array(sol[0]) 
x2=np.array(sol[1]) 
e=x1-1 
de=x2 
u=kp*e+kd*de 
for k in range(2): 
    if k==0: 
     plt.subplot(2,1,k+1) 
     plt.plot(t,sol[k],label='x1') 
     plt.plot(t,sol[k+1],label='x2') 
     plt.legend(loc='lower right') 
    else: 
     plt.subplot(2,1,k+1) 
     plt.plot(t,u) 
plt.show() 

Answer 1

ソルバーと、それが正当化されていない、訪れたタイムステップを前提としています。積分をハッキングすると、たとえそれが数学的に健全であっても（注文1の積分方法を与えるintegral = integral + e*(t-told)のように見える）、注文するだけの幸運な場合は、おそらく1までの積分方法の順序を減らしますこのシステムを実装するための2

数学的に正しい方法はすなわち、の誘導体がeである、eの積分のための第3の可変を導入することです。正しいオーダー1システムが次元3でなければならないということは、あなたのシステムに3つの微分/積分演算があるという事実を読むことができます（eを除く）。お使いのシステムが必要なグローバルな動的変数、（また、より効率的または読みやすいと思われるものは何でも、パラメータとして渡すことができます）定数のみが存在しないことを

def pid(t,Y): 
    x1, x2, x3 =Y 
    e=x1-1 
    x1dot = x2 
    edot = x1dot 
    x3dot = e 
    u=kp*e+kd*edot+ki*x3 
    x2dot=u-5*x1-2*x2 
    return[x1dot, x2dot, x3dot] 

注なっていることと。

の初期値が必要です。システムから統合変数が表示されない場合は、コードに0と表示されているようです。

Answer 2

まず、pid関数に "s"変数を含める必要があります。

：私は今見ることができる ' （S、T、Y）デフのpid ... '

Answer 3

最も簡単な解決策は、このクラスのプロパティとしてsとtoldを持つクラスを作成することです：

class PIDSolver: 
    def __init__(self) 
     self.t0=0 
     self.y0=[0.1,0.2] 
     self.kp,self.kd,self.ki=2,0.5,0.8 
     self.s,self.told=0,0 

    def pid(t,Y): 
     x1,x2=Y[0],Y[1] 
     e=x1-1 
     de=x2 
     self.s=(x1+self.s) 
     integral=self.s*(t-self.told) 
     self.told=t 
     #ie= 
     u=self.kp*e+self.kd*de+self.ki*integral 
     x1dot=x2 
     x2dot=u-5*x1-2*x2 
     return[x1dot,x2dot] 

問題の最初の部分です。ソリューションの次の部分にpidsolver = PIDSolver()を使用します。

Answer 4

set_f_params()メソッドを使用し、itz引数にリストを渡すことで、この問題を自分で解決しました。また、pid()の3番目の引数、つまりpid(t,Y,arg)を渡しました。最後に私はs,told=arg[0],arg[1]を割り当てました。