基于神經(jīng)網(wǎng)絡(luò)的預(yù)測控制模型仿真

基于神經(jīng)網(wǎng)絡(luò)的預(yù)測控制模型仿真摘要：本文利用一種權(quán)值可以在線調(diào)整的動態(tài)BP神經(jīng)網(wǎng)絡(luò)對模型預(yù)測誤差進(jìn)行擬合并與預(yù)測模型一起構(gòu)成動態(tài)組合預(yù)測器，在此根底上形成對模型誤差具有動態(tài)補償能力的預(yù)測控制算法。該算法顯著提高了預(yù)測精度，增強了預(yù)測控制算法的魯棒性。關(guān)鍵詞：預(yù)測控制神經(jīng)網(wǎng)絡(luò)動態(tài)矩陣誤差補償1.引言動態(tài)矩陣控制(DMC)是一種適用于漸近穩(wěn)定的線性或弱非線性對象的預(yù)測控制算法，目前已廣泛應(yīng)用于工業(yè)過程控制。它基于對象階躍響應(yīng)系數(shù)建立預(yù)測模型，因此建模簡單，同時采用多步滾動優(yōu)化與反應(yīng)校正相結(jié)合，能直接處理大時滯對象，并具有良好的跟蹤性能和較強的魯棒性。但是，DMC算法在實際控制中存在一系列問題，模型失配是其中普遍存在的一個問題，并會不同程度地影響系統(tǒng)性能。DMC在實際控制中產(chǎn)生模型失配的原因主要有2個，一是諸如建模誤差、環(huán)境干擾等因素，它會在實際控制的全程范圍內(nèi)引起DMC的模型失配；二是實際系統(tǒng)的非線性特性，這一特性使得被控對象的模型發(fā)生變化，此時假設(shè)用一組固定的階躍響應(yīng)數(shù)據(jù)設(shè)計控制器進(jìn)行全程范圍的控制，必然會使實際控制在對象的非建模區(qū)段內(nèi)出現(xiàn)模型失配。針對DMC模型失配問題，已有學(xué)者進(jìn)行了大量的研究，并取得了豐富的研究成果，其中有基于DMC控制參數(shù)在線辨識的智能控制算法，基于模型在線辨識的自校正控制算法以及用神經(jīng)元網(wǎng)絡(luò)進(jìn)行模型辨識、在辨識的根底上再進(jìn)行動態(tài)矩陣控制等。這些算法盡管進(jìn)行在線辨識修正對象模型參數(shù)，仍對對象降階建模誤差(結(jié)構(gòu)性建模誤差)的魯棒性不好，并對隨機噪聲干擾較敏感。針對以上問題，出現(xiàn)了基于誤差校正的動態(tài)矩陣控制算法。這些文獻(xiàn)用基于時間序列預(yù)測的數(shù)學(xué)模型誤差代替原模型誤差，得到對未來誤差的預(yù)測。有人還將這種誤差預(yù)測方法引入動態(tài)矩陣控制，并應(yīng)用于實際。這種方法雖然使系統(tǒng)表現(xiàn)出良好的穩(wěn)定性，但建立精確的誤差數(shù)學(xué)模型還存在一定的困難。本文利用神經(jīng)網(wǎng)絡(luò)通過訓(xùn)練學(xué)習(xí)能逼近任意連續(xù)有界函數(shù)的特點，建立了一種采用BP神經(jīng)網(wǎng)絡(luò)進(jìn)行預(yù)測誤差補償?shù)腄MC預(yù)測控制模型。其中神經(jīng)網(wǎng)絡(luò)預(yù)測誤差描述了在預(yù)測模型中未能包含的一切不確定性信息，可以歸結(jié)為用BP神經(jīng)網(wǎng)絡(luò)基于一系列過去的誤差信息預(yù)測未來的誤差，它作為模型預(yù)測的重要補充，不僅降低建立數(shù)學(xué)模型的負(fù)擔(dān)，而且還可以彌補在對象模型中已簡化或無法加以考慮的一切其他因素。本文通過進(jìn)行仿真，驗證了基于神經(jīng)網(wǎng)絡(luò)誤差補償?shù)念A(yù)測控制算法的有效性及優(yōu)越性，從而增強了動態(tài)矩陣控制算法的魯棒性。2.BP神經(jīng)網(wǎng)絡(luò)模型BP神經(jīng)網(wǎng)絡(luò)模型處理信息的根本原理是：輸入信號Xi通過中間節(jié)點〔隱層點〕作用于輸出節(jié)點，經(jīng)過非線形變換，產(chǎn)生輸出信號Yk，網(wǎng)絡(luò)訓(xùn)練的每個樣本包括輸入向量X和期望輸出量t，網(wǎng)絡(luò)輸出值Y與期望輸出值t之間的偏差，通過調(diào)整輸入節(jié)點與隱層節(jié)點的聯(lián)接強度取值Wij和隱層節(jié)點與輸出節(jié)點之間的聯(lián)接強度Tjk以及閾值，使誤差沿梯度方向下降，經(jīng)過反復(fù)學(xué)習(xí)訓(xùn)練，確定與最小誤差相對應(yīng)的網(wǎng)絡(luò)參數(shù)〔權(quán)值和閾值〕，訓(xùn)練即告停止。此時經(jīng)過訓(xùn)練的神經(jīng)網(wǎng)絡(luò)即能對類似樣本的輸入信息，自行處理輸出誤差最小的經(jīng)過非線形轉(zhuǎn)換的信息。2.1BP模型建立BP網(wǎng)絡(luò)模型包括其輸入輸出模型、作用函數(shù)模型、誤差計算模型和自學(xué)習(xí)模型?！?〕節(jié)點輸出模型隱節(jié)點輸出模型：Oj=f(∑Wij×Xi-qj)(1)輸出節(jié)點輸出模型：Yk=f(∑Tjk×Oj-qk)(2)f-非線形作用函數(shù)；q-神經(jīng)單元閾值?！?〕作用函數(shù)模型作用函數(shù)是反映下層輸入對上層節(jié)點刺激脈沖強度的函數(shù)，又稱刺激函數(shù)，一般取為(0,1)內(nèi)連續(xù)取值Sigmoid函數(shù)：f(x)=1/(1+e-x)〔3〕〔3〕誤差計算模型誤差計算模型是反映神經(jīng)網(wǎng)絡(luò)期望輸出與計算輸出之間誤差大小的函數(shù)：Ep=1/2×∑(tpi-Opi)2tpi-i節(jié)點的期望輸出值；Opi-i節(jié)點計算輸出值。(4)〔4〕自學(xué)習(xí)模型神經(jīng)網(wǎng)絡(luò)的學(xué)習(xí)過程，即連接下層節(jié)點和上層節(jié)點之間的權(quán)重拒陣Wij的設(shè)定和誤差修正過程。BP網(wǎng)絡(luò)有師學(xué)習(xí)方式-需要設(shè)定期望值和無師學(xué)習(xí)方式-只需輸入模式之分。自學(xué)習(xí)模型為△Wij(n+1)=h×Фi×Oj+a×△Wij(n)(5)h-學(xué)習(xí)因子；Фi-輸出節(jié)點i的計算誤差；Oj-輸出節(jié)點j的計算輸出；a-動量因子。2.2BP網(wǎng)絡(luò)模型的缺陷分析及優(yōu)化策略〔1〕學(xué)習(xí)因子h的優(yōu)化采用變步長法根據(jù)輸出誤差大小自動調(diào)整學(xué)習(xí)因子，來減少迭代次數(shù)和加快收斂速度。h=h+a×(Ep(n)-Ep(n-1))/Ep(n)a為調(diào)整步長，0~1之間取值〔6〕〔2〕隱層節(jié)點數(shù)的優(yōu)化隱層節(jié)點數(shù)的多少對網(wǎng)絡(luò)性能的影響較大，當(dāng)隱節(jié)點數(shù)太多時，會導(dǎo)致網(wǎng)絡(luò)學(xué)習(xí)時間過長，甚至不能收斂；而當(dāng)隱節(jié)點數(shù)過小時，網(wǎng)絡(luò)的容錯能力差。利用逐步回歸分析法并進(jìn)行參數(shù)的顯著性檢驗來動態(tài)刪除一些線形相關(guān)的隱節(jié)點，節(jié)點刪除標(biāo)準(zhǔn)：當(dāng)由該節(jié)點出發(fā)指向下一層節(jié)點的所有權(quán)值和閾值均落于死區(qū)〔通常取±0.1、±0.05等區(qū)間〕之中，那么該節(jié)點可刪除。最正確隱節(jié)點數(shù)L可參考下面公式計算：L=(m+n)1/2+c(7)m-輸入節(jié)點數(shù)；n-輸出節(jié)點數(shù)；c-介于1～10的常數(shù)?！?〕輸入和輸出神經(jīng)元確實定利用多元回歸分析法對神經(jīng)網(wǎng)絡(luò)的輸入?yún)?shù)進(jìn)行處理，刪除相關(guān)性強的輸入?yún)?shù)，來減少輸入節(jié)點數(shù)?！?〕算法優(yōu)化由于BP算法采用的是剃度下降法，因而易陷于局部最小并且訓(xùn)練時間較長。用基于生物免疫機制地既能全局搜索又能防止未成熟收斂的免疫遺傳算法IGA取代傳統(tǒng)BP算法來克服此缺點。2.3BP神經(jīng)網(wǎng)絡(luò)仿真假設(shè)備控對象的近似數(shù)學(xué)模型為：〔8〕式中，系統(tǒng)a〔k〕是慢時變的，。神經(jīng)網(wǎng)絡(luò)的結(jié)構(gòu)選擇2-6-1，學(xué)習(xí)速度η=0.5和慣性系數(shù)α=0.05，輸入指令信號為：rin(k)=0.5*sin(6*л*t)。仿真程序：clearall;closeall;xite=0.50;alfa=0.05;w2=rand(6,1);w2_1=w2;w2_2=w2;w1=rand(2,6);w1_1=w1;w1_2=w1;dw1=0*w1;x=[0,0]';u_1=0;y_1=0;I=[0,0,0,0,0,0]';Iout=[0,0,0,0,0,0]';FI=[0,0,0,0,0,0]';ts=0.001;fork=1:1:1000time(k)=k*ts;u(k)=0.50*sin(3*2*pi*k*ts);a(k)=1.2*(1-0.8*exp(-0.1*k));y(k)=a(k)*y_1/(1+y_1^2)+u_1;forj=1:1:6I(j)=x'*w1(:,j);Iout(j)=1/(1+exp(-I(j)));endyn(k)=w2'*Iout;e(k)=y(k)-yn(k);w2=w2_1+(xite*e(k))*Iout+alfa*(w2_1-w2_2);forj=1:1:6FI(j)=exp(-I(j))/(1+exp(-I(j)))^2;endfori=1:1:2forj=1:1:6dw1(i,j)=e(k)*xite*FI(j)*w2(j)*x(i);endendw1=w1_1+dw1+alfa*(w1_1-w1_2);x(1)=u(k);x(2)=y(k);w1_2=w1_1;w1_1=w1;w2_2=w2_1;w2_1=w2;u_1=u(k);y_1=y(k);endsubplot(2,1,1);plot(time,y,'r',time,yn,'b');xlabel('times');ylabel('實際輸出和仿真輸出');gridonsubplot(2,1,2);plot(time,y-yn,'r');xlabel('times');ylabel('error');gridon圖2-1神經(jīng)網(wǎng)絡(luò)仿真圖3．動態(tài)矩陣控制DMC算法分為三局部:預(yù)測模型、滾動優(yōu)化和反應(yīng)校正。3.1預(yù)測模型在DMC算法中，首先需要測定對象單位階躍響應(yīng)的采樣值,i=1,2,……。其中，T為采樣周期。對于漸進(jìn)穩(wěn)定的對象，階躍響應(yīng)在某一，后將趨于平穩(wěn)，以至(i>N)與的誤差和量化誤差及測量誤差有相同的數(shù)量級，因而可認(rèn)為，已近似等于階躍響應(yīng)的穩(wěn)態(tài)值。這樣，對象的動態(tài)信息就可以近似用有限集合加以描述。這個集合的參數(shù)構(gòu)成了DMC的模型參數(shù)，向量a=稱為模型向量，N稱為模型時域長度。雖然階躍響應(yīng)是一種非參數(shù)模型，但由于線性系統(tǒng)具有比例和疊加性質(zhì)，故利用這組模型參數(shù)，已足以預(yù)測在任意輸入作用下系統(tǒng)在未來時刻的輸出值。在t=kT時刻，假設(shè)控制量不再變化時系統(tǒng)在未來N個時刻的輸出值為,那么，在控制增量作用后系統(tǒng)的輸出可由〔3-1〕預(yù)測，其中表示在t=kT時刻預(yù)測的尚無作用時未來N個時刻的系統(tǒng)輸出。表示在t=kT時刻預(yù)測的有控制增量作用時未來N個時刻的系統(tǒng)輸出。為階躍響應(yīng)模型向量，其元素為描述系統(tǒng)動態(tài)特性的N個階躍響應(yīng)系數(shù)。式中，上標(biāo)~表示預(yù)測，表示在t=kT時刻預(yù)測t=(k+i)T時刻。同樣，如果考慮到現(xiàn)在和未來M個時刻控制增量的變化，在t=kT時刻預(yù)測在控制增量,?,作用下系統(tǒng)在未來P個時刻的輸出為〔3-2〕式中為t=kT時刻預(yù)測的無控制增量時未來P個時刻的系統(tǒng)輸出。為t=kT時刻預(yù)測的有M個控制增量,……，時未來P個時刻的系統(tǒng)輸出。為從現(xiàn)在起M個時刻的控制增量。稱為動態(tài)矩陣，其元素為描述系統(tǒng)動態(tài)特性的階躍響應(yīng)系數(shù)。3.2滾動優(yōu)化DMC是一種以優(yōu)化確定控制策略的算法。在采樣時刻t=kT的優(yōu)化性能指標(biāo)可取為+〔3-3〕即通過選擇該時刻起M個時刻的控制增量,?，，使系統(tǒng)在未來P〔〕個時刻的輸出值，…，盡可能接近其期望值，…。性能指標(biāo)中的第二項是對控制增量的約束，即不允許控制量的變化過于劇烈。式中，，為權(quán)系數(shù)，P和M分別稱為優(yōu)化時域長度和控制時域長度。顯然，在不同時刻，優(yōu)化性能指標(biāo)是不同的，但其相對形式卻是一致的，都具有類似于(3-1)的形式，所謂“滾動優(yōu)化〞，就是指優(yōu)化時域隨時間不斷地向前推移。引入向量和矩陣記號，Q=diag(,…，)，R=diag(,…，)那么優(yōu)化性能指標(biāo)式(2-1)可改寫為(3-4)式中，Q,R分別稱為誤差權(quán)矩陣和控制權(quán)矩陣。在不考慮輸入輸出約束的情況下，在t=kT時刻，，均為己知，使J(k)取最小的可通過極值必要條件求得(3-5)這就是t=kT時刻解得的最優(yōu)控制增量序列。由于這一最優(yōu)解完全是基于預(yù)測模型求得的因而是開環(huán)最優(yōu)解.3.3反應(yīng)校正由于模型誤差、弱非線性特性及其它在實際過程中存在的不確定因素，按預(yù)測模型式(3-2)得到的開環(huán)最優(yōu)控制規(guī)律式(3-5)不一定能導(dǎo)致系統(tǒng)輸出緊密地跟隨期望值，它也不能顧及對象受到的擾動。為了糾正模型預(yù)測與實際的不一致，必須及時地利用過程的誤差信息對輸出預(yù)測值進(jìn)行修正，而不應(yīng)等到這M個控制增量都實施后再作校正。為此，在t=kT時刻首先實施中的第一個控制作用=(3-6)(3-7)其中(3-8)由于已作用于對象，對系統(tǒng)未來輸出的預(yù)測便要疊加上產(chǎn)生的影響，即由式(2-1)算出。到下一個采樣時刻t=(k+1)T，不是繼續(xù)實施最優(yōu)解中的第二個分量，而是檢測系統(tǒng)的實際輸出y(k+1)，并與按模型預(yù)測算得的該時刻輸出，即中的第一個分量進(jìn)行比擬，構(gòu)成預(yù)測誤差。(3-9)這一誤差反映了模型中未包含的各種不確定因素，如模型失配、干擾等。由于預(yù)測誤差的存在，以后各時刻輸出值的預(yù)側(cè)也應(yīng)在模型預(yù)測的根底上加以校正，這些未來誤差的預(yù)測，可通過對現(xiàn)時誤差e(k+1)加權(quán)系數(shù)h,(i=1,2,…,N)得到〔3-10〕式中為t=(k+1)T時刻經(jīng)誤差校正后所預(yù)測的系統(tǒng)在t=(k+i)T〔i=1,…，N〕時刻的輸出.為誤差校正向量，其中，=1。經(jīng)校正后的的各分量中。除第一項外.其余各項分別是t=(k+1)T時刻在尚無等未來控制增量作用時對輸出在t=(k+2)T,…,(k+N)T時刻的預(yù)測值，它們可作為t=(k+1)T時刻的前N一1個分量，即,i=1,…，N-1而中的最后一個分量。即t=(k+1)T時刻對i=(k+1+N)T輸出的預(yù)測，可由來近似，即=，上述關(guān)系可用向量形式表示(3-11)其中為移位矩陣。在t=(k+1)T時刻，有了,就又可以像上面所述t=kT時刻那樣進(jìn)行新的預(yù)測優(yōu)化，整個控制就是在這樣推移的過程中滾動進(jìn)行.由此可以看到，整個動態(tài)矩陣控制算法是由預(yù)測、控制、校正三局部組成的，該算法結(jié)構(gòu)可用圖2-1加以描述。圖中粗箭頭表示向量流，細(xì)箭頭表示純量流。在每一個采樣時刻未來P個時刻的期望輸出與預(yù)測輸出所構(gòu)成的偏差向量按式(3-6)與動態(tài)向量點乘。得到該時刻的控制增量。這一控制增量一方面通過數(shù)字積分(累加)運算求出控制量u(k)作用于對象;另一方面與階躍響應(yīng)向量a相乘，并按式(3-1)計算出在其作用后所預(yù)測的系統(tǒng)輸出。到了下一個采樣時刻，首先測定系統(tǒng)的實際輸出y(k+1)，并與原來預(yù)測的該時刻的值相比擬。按(3-9)算出預(yù)測誤差e(k+1)。這一誤差與校正向量h相乘后，再按式(3-10)校正預(yù)測的輸出值。由于時間的推移，經(jīng)校正的預(yù)測輸出將按式(3-11)移位，并置定為該時刻的預(yù)測初值.圖中的，表示時移算子，如果把新的時刻重新定義為k時刻，那么預(yù)測初值的前P個分量將與期望輸出一起，參與新時刻控制增量的計算。如此循環(huán)，整個過程將反復(fù)在線進(jìn)行。DMC算法是一種增量算法?？梢宰C明，不管有否模型誤差，它總能將系統(tǒng)輸出調(diào)節(jié)到期望值而不產(chǎn)生靜差。對于作用在對象輸入端的階躍形式的擾動，該算法也總能使系統(tǒng)輸出回復(fù)到原來的設(shè)定狀態(tài).3.4有約束多變量動態(tài)矩陣控制及其線性化不失一般性,假設(shè)階躍響應(yīng)的截斷時域長度均為N,控制時域長度均為M,預(yù)測時域長度均為P?；陔A躍響應(yīng)的預(yù)測模型為(3-12)其中為k時刻第j個輸入對第i個輸出的奉獻(xiàn)，為第i個輸出對第j個輸入的第s個階躍響應(yīng)系數(shù)。那么第i個模型輸出為〔3-13〕令，考慮到反應(yīng)校正作用，把控制量分解為當(dāng)前及將來未知局部和過去局部。可得k時刻輸出的P步預(yù)測為：〔3-14〕其中，和分別為和階矩陣，其元素由決定，為維校正列向量，而將化成全量形式，有〔3-15〕，將式〔3-15〕代入〔3-14〕得〔3-16〕令，那么式〔2-16〕成為〔3-17〕設(shè)第i個過程的P步參數(shù)輸出為為跟蹤參考軌跡，應(yīng)使為使上式性能指標(biāo)線性化，將其改寫成〔3-18〕其中，和表示相應(yīng)維數(shù)的正、負(fù)偏差變量列向量，對此二向量求其最小指對其范數(shù)求最小。令結(jié)合式(2-18)和〔2-17〕可得：〔3-19〕對控制施加幅值約束，不失一般性，假設(shè)控制幅值的下限為零，即〔3-20〕其中，為第j個控制的最大幅值，。如果控制幅值的下限約束為負(fù)，那么可通過坐標(biāo)或利用線性規(guī)劃變量非負(fù)化的方法使其為零，以滿足單純形法迭代尋優(yōu)的要求。對式〔3-20〕引入相應(yīng)維數(shù)的正、負(fù)偏差變量列向量，使不希望偏差最小化，那么式〔3-20〕為：〔3-21〕對控制變化速率施加上、下限約束,即〔3-22〕其中，，為第j個控制變化速率的上、下限。同理,引入相應(yīng)維數(shù)的正、負(fù)偏差變量列向量,最小化不希望偏差,那么式(2-22)可表示為〔3-23〕3.5動態(tài)矩陣控制仿真設(shè)備控對象的近似數(shù)學(xué)模型為：式中，系統(tǒng)a〔k〕是慢時變的，，輸入指令信號為：rin(k)=10。仿真程序：clearall;

closeall;u_1=0;y_1=0;fork=1:1:1000u(k)=1;a(k)=1.2*(1-0.8*exp(-0.1*k));yg(k)=a(k)*y_1/(1+y_1^2)+u_1;y0(k)=yg(k);y_1=yg(k);u_1=u(k);endTs=5;

N=200;P=25;m=20;Sv=10;yg_1=0;A=zeros(P,m);a=zeros(N,1);

fori=1:N

a(i)=y0(i);

end

fori=1:P

forj=1:m

ifi-j+1>0

A(i,j)=a(i-j+1);

end

end

end

K=inv(A'*A+eye(m))*A';y=zeros(N,1);

u=zeros(N,1);

e=zeros(N,1);A0=zeros(P,N-1);

fori=1:P

forj=N-2:-1:1

ifN-j+1+i-1<=N

A0(i,j)=a(N-j+1+i-1)-a(N-j+i-1);

else

A0(i,j)=0;

end

end

A0(i,N-1)=a(i+1);

end

fork=2:N

Uk_1=zeros(N-1,1);

fori=1:N-1

ifk-N+i<=0

Uk_1(i)=0;

else

Uk_1(i)=u(k-N+i);

end

end

Y0=A0*Uk_1;a(k)=1.2*(1-0.8*exp(-0.1*k));yg(k)=a(k)*y_1/(1+y_1^2)+u_1;u_1=u(k);yg_1=yg(k);

e(k)=y(k)-yg(k);

Ysk=zeros(P,1);

fori=1:P

Ysk(i)=Sv;

end

Ek=zeros(P,1);

fori=1:P

Ek(i)=0.1*i*e(k);

end

dertu=K*(Ysk-Y0-Ek);

fori=1:m

ifk+i-1<=N

u(k+i-1)=u(k+i-1-1)+dertu(i);

end

end

temp=0;

forj=1:N-1

ifk-j<=0

temp;

else

ifk-j-1<=0

temp=temp+a(j)*u(k-j);

else

temp=temp+a(j)*(u(k-j)-u(k-j-1));

end

end

end

ifk-N<=0

y(k)=temp+e(N);

else

y(k)=temp+a(N)*u(k-N)+e(N);

end

end

t=Ts.*(1:N);

subplot(2,1,1);

plot(t,y,'r');

legend('輸出曲線y');

title('輸出曲線');

xlabel('Time')

ylabel('振幅')

gridon

subplot(2,1,2);

plot(t,u,'g');

legend('控制作用u')

title('控制作用');

xlabel('Time')

ylabel('振幅')

gridon

clc;圖3-1動態(tài)矩陣控制效果圖4.BP網(wǎng)絡(luò)誤差補償預(yù)測模型仿真DMC是基于一個線性模型對過程的輸出進(jìn)行預(yù)測的，雖然對模型失配有較強的魯棒性，但對隨機突發(fā)干擾以及實際系統(tǒng)中的非線性，仍然存在著難以及時控制的缺乏，從而產(chǎn)生預(yù)測偏差，所以采用BP神經(jīng)網(wǎng)絡(luò)計算預(yù)測偏差，從而進(jìn)行補償假設(shè)備控對象的近似數(shù)學(xué)模型為：式中，系統(tǒng)a〔k〕是慢時變的，。選用第二章的BP神經(jīng)網(wǎng)絡(luò)，學(xué)習(xí)速度η=0.5和慣性系數(shù)α=0.5，輸入指令信號為：rin(k)=10。仿真程序：clearall;closeall;xite=0.50;alfa=0.05;w2=rand(6,1);w2_1=w2;w2_2=w2;w1=rand(2,6);w1_1=w1;w1_2=w1;dw1=0*w1;x=[0,0]';u_1=0;y_1=0;I=[0,0,0,0,0,0]';Iout=[0,0,0,0,0,0]';FI=[0,0,0,0,0,0]';ts=0.001;fork=1:1:1000time(k)=k*ts;u(k)=0.50*sin(3*2*pi*k*ts);a(k)=1.2*(1-0.8*exp(-0.1*k));y(k)=a(k)*y_1/(1+y_1^2)+u_1;forj=1:1:6I(j)=x'*w1(:,j);Iout(j)=1/(1+exp(-I(j)));endyn(k)=w2'*Iout;e(k)=y(k)-yn(k);w2=w2_1+(xite*e(k))*Iout+alfa*(w2_1-w2_2);forj=1:1:6FI(j)=exp(-I(j))/(1+exp(-I(j)))^2;endfori=1:1:2forj=1:1:6dw1(i,j)=e(k)*xite*FI(j)*w2(j)*x(i);endendw1=w1_1+dw1+alfa*(w1_1-w1_2);x(1)=u(k);x(2)=y(k);w1_2=w1_1;w1_1=w1;w2_2=w2_1;w2_1=w2;u_1=u(k);y_1=y(k);endfigure(1);plot(time,y,'r',time,yn,'b');xlabel('times');ylabel('yandyn');gridonts=0.1;fork=1:1:200time(k)=k*ts;u(k)=1;a(k)=1.2*(1-0.8*exp(-0.1*k));y(k)=a(k)*y_1/(1+y_1^2)+u_1;forj=1:1:6I(j)=x'*w1(:,j);Iout(j)=1/(1+exp(-I(j)));endyn(k)=w2'*Iout;e(k)=y(k)-yn(k);w2=w2_1+(xite*e(k))*Iout+alfa*(w2_1-w2_2);forj=1:1:6FI(j)=exp(-I(j))/(1+exp(-I(j)))^2;endfori=1:1:2forj=1:1:6dw1(i,j)=e(k)*xite*FI(j)*w2(j)*x(i);endendw1=w1_1+dw1+alfa*(w1_1-w1_2);x(1)=u(k);x(2)=y(k);w1_2=w1_1;w1_1=w1;w2_2=w2_1;w2_1=w2;u_1=u(k);y_1=y(k);y0(k)=yn(k);endTs=5;N=200;P=25;m=20;Sv=10;yg_1=0;A=zeros(P,m);a=zeros(N,1);fori=1:Na(i)=y0(i);endfori=1:Pforj=1:mifi-j+1>0A(i,j)=a(i-j+1);endendendK=inv(A'*A+eye(m))*A';y=zeros(N,1);u=zeros(N,1);e=zeros(N,1);A0=zeros(P,N-1);fori=1:Pforj=N-2:-1:1ifN-j+1+i-1<=NA0(i,j)=a(N-j+1+i-1)-a(N-j+i-1);elseA0(i,j)=0;endendA0(i,N-1)=a(i+1);endfork=2:NUk_1=zeros(N-1,1);fori=1:N-1ifk-N+i<=0Uk_1(i)=0;elseUk_1(i)=u(k-N+i);endendY0=A0*Uk_1;a(k)=1.2*(1-0.8*exp(-0.1*k));yg(k)=a(k)*y_1/(1+y_1^2)+u_1;u_1=u(k);yg_1=yg(k);e(k)=y(k)-yg(k);Ysk=zeros(P,1);fori=1:PYsk(i)=Sv;endEk=zeros(P,1);fori=1:PEk(i)=0.1*i*e(k);enddertu=K*(Ysk-Y0-Ek);fori=1:mifk+i-1<=Nu(k+i-1)=u(k+i-1-1)+dertu(i);endendtemp=0;forj=1:N-1ifk-j<=0temp;elseifk-j-1<=0temp=temp+a(j)*u(k-j);elsetemp=temp+a(j)*(u(k-j)-u(k-j-1));endendendifk-N<=0y(k)=temp+e(N);elsey(k)=temp+a(N)*u(k-N)+e(N);endendt=Ts.*(1:N);figure(3);subplot(2,1,1);plot(t,y,'b',t,yg,'r');legend('預(yù)測輸出曲線');title('輸出曲線');xlabel('Time')ylabel('振幅')gridonsubplot(2,1,2);plot(t,u,'g');legend('控制作用u')title('控制作用');xlabel('Time')ylabel('振幅')gridonxite=0.50;alfa=0.05;w2=rand(10,1);w2_1=w2;w2_2=w2;w1=rand(6,10);w1_1=w1;w1_2=w1;dw1=0*w1;u_1=0;y_1=0;I=zeros(10,1);Iout=zeros(10,1);FI=zeros(10,1);x=[0,0,0,0,0,0]';fork=3:Nx(1)=yg(k);x(2)=yg(k-1);x(3)=yg(k-2);x(4)=u(k);x(5)=u(k-1);x(6)=u(k-2);forj=1:1:10I(j)=x'*w1(:,j);Iout(j)=1/(1+exp(-I(j)));endes(k)=w2'*Iout;e1(k)=es(k)-e(k);w2=w2_1-(xite*0.1*e1(k))*Iout+alfa*(w2_1-w2_2);forj=1:1:10FI(j)=exp(-I(j))/(1+exp(-I(j)))^2;endfori=1:1:6forj=1:1:10dw1(i,j)=e1(k)*xite*FI(j)*w2(j)*x(i);endendw1=w1_1+dw1+alfa*(w1_1-w1_2);w1_2=w1_1;w1_1=w1;w2_2=w2_1;w2_1=w2;endt=Ts.*(1:N);figure(2);plot(t,es,'r',t,e,'b');xlabel('times');ylabel('esande');gridonfork=3:NUk_1=zeros(N-1,1);fori=1:N-1ifk-N+i<=0Uk_1(i)=0;