基于Matlab仿真的線性與非線性動態(tài)系統(tǒng)建模教學方法探討

王直杰

[摘要]本文采用Matlab仿真產(chǎn)生模擬的可用于線性與非線性動態(tài)系統(tǒng)建模的輸入輸出數(shù)據(jù)。根據(jù)這些數(shù)據(jù)，設計程序生成建模所需的數(shù)據(jù)矩陣。同時給出了最小二乘法、最大似然估計、BP神經(jīng)網(wǎng)絡、RBF神經(jīng)網(wǎng)絡建模的樣例程序。為學生熟悉、比較及應用各種線性與非線性動態(tài)系統(tǒng)建模方法提供了條件。

[關鍵詞]動態(tài)系統(tǒng)建模 仿真 人工神經(jīng)網(wǎng)絡

在講授完線性與非線性動態(tài)系統(tǒng)建模方法基本理論以后，需要讓學生進行上機實驗[1]以達到以下目的：1進一步加深理解學習的基本理論；2各種參數(shù)的選擇對建模效果的影響；3各種建模方法的優(yōu)缺點及適用場合。為了達到這些目的，我們模擬實際應用時的情況，提供給學生模擬的輸入輸出數(shù)據(jù)，以便學生利用這些數(shù)據(jù)，進行編程建立模型。同時我們編程實現(xiàn)基于最小二乘法、最大似然估計、BP神經(jīng)網(wǎng)絡、RBF神經(jīng)網(wǎng)絡的動態(tài)系統(tǒng)建模方法[2，3]，學生可以利用這些Matlab程序進行各種方法的學習、各種情況下建模效果的對比，以及各種方法的適用場合的對比。

一、 基于Matlab仿真的線性動態(tài)系統(tǒng)最小二乘法建模的教學

假設系統(tǒng)的差分方程為：y（k）=-a1y（k-1）-a2y（k-2）-…-any（k-n）+b0u（k）+b1u（k-1）+…+bmu（k-m）+e（k）。其中y（k）為輸出，u（k）為輸入，e（k）為模型殘差。假定建模用的數(shù)據(jù)序列從y（k）開始，則構建以下數(shù)據(jù)矩陣及數(shù)據(jù)向量：

如果模型殘差為白噪聲（實際情況多為有色噪聲，但當噪聲強度不大時，可近似當作白噪聲處理），則根據(jù)最小二乘法，由這 組數(shù)據(jù)估計得到的參數(shù) 。

以上的最小二乘法需要輸入（u（K-m），u（K-1），…，u（K+N-1））和輸出（y（K-n），y（K-n+1），…，y（K+N-1））數(shù)據(jù)，我們可用以下的Matlab程序（程序1）產(chǎn)生模擬的輸入輸出數(shù)據(jù)（程序中采用了一個簡單的二階離散系統(tǒng)，學生實驗時可換成需要的模型），并形成數(shù)據(jù)矩陣及數(shù)據(jù)向量。

程序1：

clear all

K=3；%殘差序列開始序號

N=1000；%共N組數(shù)據(jù)

n=2；%對應

m=2；%對應

u=randn（1，K+N-1）*0.1；%隨機產(chǎn)生輸入數(shù)據(jù)

for i=1：1：n

y（i）=0；%設定初始狀態(tài)

end

for k=n+1：1：K+N-1

y（k）=0.2*y（k-1）+0.5*y（k-2）+u（k）+0.3*u（k-1）；%根據(jù)輸入計算輸出

end

y=y+randn（1，K+N-1）*0.001； %在求得的單位階躍響應上疊加噪聲用以模擬測量誤差

X=zeros（N，n+m+1）；

for i=K：1：K+N-1

for j=1：1：n

X（i-K+1，j）=y（i-j）；%構成矩陣X

end

for j=n+1：1：n+m+1

X（i-K+1，j）=u（i-j+n+1）； %構成矩陣X

end

end

for i=K：1：K+N-1

Y（i-K+1）=y（i）；%構成向量Y

end

save X X；

save Y Y；

save N N；

運行程序1將生產(chǎn)數(shù)據(jù)矩陣X（保存在X.mat中）及數(shù)據(jù)向量Y（保存在Y.mat中）。以下是根據(jù)最小二乘法估計參數(shù)的程序（程序2）。

程序2：

clear all

load X；

load Y；

XT=X'；

sita=（XT*X）^（-1）*XT*Y'%計算得到估計的參數(shù)sita

運行程序2，將得到估計的參數(shù)sita。如某一次運行中sita=[0.2004，0.4998，0.9998，0.2995，-0.0005]T，與模型中的值（見程序1）a1=0.2，a2=0.5，b0=1，b1=0.3，b2=0非常接近。

二、 基于Matlab仿真的線性動態(tài)系統(tǒng)最大似然估計建模的教學

在模型殘差為白噪聲的假設下，最大似然估計和最小二乘法在估計參數(shù)時是相同的，但最大似然估計還能估計出噪聲的強度，程序3為相應的程序。

程序3：

clear all

load X；

load Y；

load N；

XT=X'；

sita=（XT*X）^（-1）*XT*Y'

Z=Y'-X*sita；

StdV=sqrt（Z'*Z/N）%估計標準差

運行程序3，得到和程序2一樣的參數(shù)估計值，除此之外，還能得到噪聲的標準差（或方差）。如某一次的運行結果為StdV=0.0011，和模型中值（見程序1）0.001非常接近。

三、 基于Matlab仿真的非線性動態(tài)系統(tǒng)BP神經(jīng)網(wǎng)絡建模的教學

基于神經(jīng)網(wǎng)絡的非線性動態(tài)系統(tǒng)建模的輸入輸出的樣本數(shù)據(jù)組織如圖1所示。

圖1：基于神經(jīng)網(wǎng)絡的動態(tài)系統(tǒng)建模的樣本數(shù)據(jù)的組織

如圖1所示，當輸入為y（K-1），y（K-2），…，y（K-n），u（K），u（K-1），…，u（K-m）時，期望輸出為y（K），因此BP神經(jīng)網(wǎng)絡的輸入輸出數(shù)據(jù)樣本對為（x，y），其中x，y分別為上述最小二乘法中的數(shù)據(jù)矩陣及數(shù)據(jù)向量，將程序1中的動態(tài)系統(tǒng)以非線性動態(tài)系統(tǒng)代替（如將y（k）=0.2*y（k-1）+0.5*y（k-2）+u（k）+0.3*u（k-1）

改為y（k）=0.2*y（k-1）+0.5*y（k-2）+2*u（k）*u（k）+0.3*u（k-1）），采用程序1產(chǎn)生非線性動態(tài)系統(tǒng)的數(shù)據(jù)，然后設計如下的基于BP神經(jīng)網(wǎng)絡的非線性動態(tài)系統(tǒng)建模程序（程序4）。

程序4：

clear all

load X；

load Y；

net = newff（X'，Y，10）；

net.trainParam.epochs = 100；

net.trainParam.goal = 0.000001；

net = train（net，X'，Y）；

Y1 = sim（net，X'）；

plot（Y，'s-'）；

hold on

plot（Y1，'*-'）；

save net net；

從程序4運行后的產(chǎn)生的圖形中可以看出BP神經(jīng)網(wǎng)絡訓練的效果。

四、 基于RBF神經(jīng)網(wǎng)絡的非線性動態(tài)系統(tǒng)建模教學

基于RBF神經(jīng)網(wǎng)絡的非線性動態(tài)系統(tǒng)建模的樣本數(shù)據(jù)的組織和BP神經(jīng)網(wǎng)絡相同，因此可以使用BP神經(jīng)網(wǎng)絡建模時所用的樣本數(shù)據(jù)進行建模實驗，以便對比兩種網(wǎng)絡的建模效果。以下（程序5）是RBF神經(jīng)網(wǎng)絡建模的樣例程序。

程序5：

clear all

load X；

load Y；

net=newrb（X'，Y，0.000001）；

Y1=sim（net，X'）；

plot（Y，'s-'）；

hold on

plot（Y1，'*-'）；

save net net；

從程序5運行后的產(chǎn)生的圖形中可以看出RBF神經(jīng)網(wǎng)絡訓練的效果。

五、 結論

本文設計了matlab程序，模擬產(chǎn)生較逼真的輸入輸出數(shù)據(jù)樣本數(shù)據(jù)供學生使用，學生可以利用這些數(shù)據(jù)進行線性與非線性動態(tài)系統(tǒng)建模實驗。同時也提供了面向動態(tài)系統(tǒng)建模的最小二乘法、最大似然估計、BP神經(jīng)網(wǎng)絡、RBF神經(jīng)網(wǎng)絡樣例程序供學生學習使用。學生可以利用這些Matlab程序進行各種方法的學習、各種情況下建模效果的對比，以及各種方法的適用場合的對比。學生也可以參照這些程序編制更加復雜的程序以解決實際的系統(tǒng)建模問題。

基金資助：本文系東華大學信息學院教改項目的研究成果。

[參考文獻]

[1]劉娣許有熊林健，基于MATLAB的“系統(tǒng)辨識”課程實驗教學改革[J].中國電力教育，2013（1）：139-140.

[2]王秀峰，盧桂章.系統(tǒng)建模與辨識[M].北京：電子工業(yè)出版社，2004.

[3]陳宗海.過程系統(tǒng)建模與仿真[M].合肥：中國科技大學出版社，1997.

（作者單位：東華大學 信息科學與技術學院 上海）