搬磚的攻城獅，愛(ài)好系統(tǒng)仿真設(shè)計(jì)
ID:0000401、系統(tǒng)的思想
系統(tǒng)的思想從古到今一直發(fā)生著演變。中國(guó)古代文化中的系統(tǒng)思想主要體現(xiàn)在對(duì)世界的認(rèn)識(shí)和解釋上，如《道德經(jīng)》中的“道生一，一生二，二生三，三生萬(wàn)物”，論證了世間萬(wàn)物由要素組成，要素之間存在相互聯(lián)系相互作用；又如春秋戰(zhàn)國(guó)諸子百家對(duì)軍事、自然、社會(huì)等現(xiàn)象進(jìn)行討論，形成了中國(guó)獨(dú)有的軍事、自然科學(xué)和社會(huì)倫理體系，并將其應(yīng)用于具體的工程項(xiàng)目。近現(xiàn)代西方的系統(tǒng)思想則更關(guān)注系統(tǒng)的運(yùn)行規(guī)律，其數(shù)學(xué)、物理、化學(xué)、生物學(xué)等基礎(chǔ)科學(xué)和自然科學(xué)的發(fā)展，使人類(lèi)對(duì)自然系統(tǒng)的組成、要素關(guān)系和運(yùn)行規(guī)律等方面的認(rèn)識(shí)取得了突破性成果，發(fā)展和完善了科學(xué)的方法論。
2、系統(tǒng)的定義
這里引用錢(qián)學(xué)森的觀點(diǎn)給出“系統(tǒng)”的定義：系統(tǒng)是由相互作用相互依賴(lài)的若干部分結(jié)合而成的、具有特定功能的有機(jī)整體，而且這個(gè)有機(jī)整體又是它從屬的更大系統(tǒng)的組成部分。由此可以提煉出認(rèn)識(shí)一個(gè)系統(tǒng)的關(guān)鍵點(diǎn)在于組成系統(tǒng)的要素和各要素如何相互作用。根據(jù)系統(tǒng)層次結(jié)構(gòu)的復(fù)雜程度，可以將系統(tǒng)分為簡(jiǎn)單系統(tǒng)和復(fù)雜系統(tǒng)：直接由要素構(gòu)成的兩層次系統(tǒng)稱(chēng)為簡(jiǎn)單系統(tǒng)，由多層子系統(tǒng)按照一定結(jié)構(gòu)構(gòu)成的系統(tǒng)稱(chēng)為復(fù)雜系統(tǒng)。例如，中國(guó)載人航天工程就是一個(gè)復(fù)雜系統(tǒng)，由航天員系統(tǒng)、空間應(yīng)用系統(tǒng)、載人飛船系統(tǒng)、貨運(yùn)飛船系統(tǒng)、長(zhǎng)征二號(hào)F運(yùn)載火箭系統(tǒng)、長(zhǎng)征七號(hào)運(yùn)載火箭系統(tǒng)、長(zhǎng)征五號(hào)B運(yùn)載火箭系統(tǒng)、酒泉發(fā)射場(chǎng)系統(tǒng)、文昌發(fā)射場(chǎng)系統(tǒng)、測(cè)控通信系統(tǒng)、空間實(shí)驗(yàn)室系統(tǒng)、空間站系統(tǒng)、著陸場(chǎng)系統(tǒng)和光學(xué)艙系統(tǒng)等十四大系統(tǒng)組成，各個(gè)系統(tǒng)又由若干個(gè)子系統(tǒng)構(gòu)成，如測(cè)控通信系統(tǒng)就包含雷達(dá)系統(tǒng)、光學(xué)系統(tǒng)等子系統(tǒng)。
圖1 中國(guó)載人航天工程系統(tǒng)組成
此外，按照時(shí)間和狀態(tài)變量的取值特性，可以將系統(tǒng)分為連續(xù)、離散系統(tǒng)和混合系統(tǒng)狀態(tài)變量隨時(shí)間連續(xù)變化的系統(tǒng)稱(chēng)為連續(xù)系統(tǒng)；時(shí)間變量和狀態(tài)變量均只取有限個(gè)離散值的系統(tǒng)稱(chēng)為離散系統(tǒng)；兼具連續(xù)和離散成分的系統(tǒng)稱(chēng)為混合系統(tǒng)。按照系統(tǒng)要素之間的關(guān)系或模型的形式，可以將系統(tǒng)分為線(xiàn)性系統(tǒng)和非線(xiàn)性系統(tǒng)：要素之間的關(guān)系能夠用線(xiàn)性數(shù)學(xué)模型描述的系統(tǒng)稱(chēng)為線(xiàn)性系統(tǒng)，否則稱(chēng)為非線(xiàn)性系統(tǒng)。除此之外，還有其他的對(duì)系統(tǒng)的分類(lèi)方式，如封閉系統(tǒng)和開(kāi)放系統(tǒng)、動(dòng)態(tài)系統(tǒng)和靜態(tài)系統(tǒng)等。
[1]中國(guó)載人航天官方網(wǎng)站 中國(guó)載人航天工程網(wǎng) http://www.cmse.gov.cn/gygc/xtzc/htyxt/。
[2] 錢(qián)學(xué)森．論宏觀建筑與與微觀建筑：杭州出版社，2001
為了研究、分析、設(shè)計(jì)和實(shí)現(xiàn)一個(gè)系統(tǒng)，必然需要進(jìn)行對(duì)系統(tǒng)的試驗(yàn)。試驗(yàn)的方法可以分為兩大類(lèi)：一種是直接在真實(shí)系統(tǒng)上進(jìn)行試驗(yàn)，一種則是根據(jù)系統(tǒng)的組成要素和要素間的相互作用，將系統(tǒng)抽象為模型，對(duì)模型進(jìn)行試驗(yàn)。系統(tǒng)模型通常分為物理模型、概念模型和數(shù)學(xué)模型。本文重點(diǎn)論述基于數(shù)學(xué)模型描述的系統(tǒng)。
通過(guò)對(duì)系統(tǒng)建模進(jìn)行試驗(yàn)和結(jié)果分析，就是所謂的“系統(tǒng)仿真”。直接在真實(shí)系統(tǒng)上進(jìn)行試驗(yàn)，顯然更有助于系統(tǒng)的研究分析，而系統(tǒng)仿真不可能完全符合系統(tǒng)實(shí)際運(yùn)行情況，那么為什么要進(jìn)行系統(tǒng)仿真呢？一言以蔽之，系統(tǒng)仿真為系統(tǒng)設(shè)計(jì)和系統(tǒng)分析提供了成本較低的有力支持。
以載人航天工程為例，其一，在系統(tǒng)設(shè)計(jì)階段，真實(shí)的系統(tǒng)尚未建立，如何了解和驗(yàn)證系統(tǒng)的性能？載人航天作為舉國(guó)體制下的重要工程，牽一發(fā)而動(dòng)全身，在給出最終系統(tǒng)方案之前，必然要先確定所采用方案中系統(tǒng)的性能，但若為每一個(gè)方案都搭建一套真實(shí)的系統(tǒng)，在資源和人力上顯然是不可承受的，系統(tǒng)仿真在這種情況下就為系統(tǒng)設(shè)計(jì)提供了有力支持；其二，在真實(shí)系統(tǒng)上進(jìn)行試驗(yàn)可能引起系統(tǒng)破壞或故障，如果用真實(shí)的飛船搭載火箭進(jìn)行發(fā)射試驗(yàn)，過(guò)程中不可避免會(huì)對(duì)系統(tǒng)產(chǎn)生破壞，其成本必然是高昂的，通過(guò)系統(tǒng)仿真可以減少試驗(yàn)成本；其三，在實(shí)際的發(fā)射過(guò)程會(huì)遇到各種復(fù)雜的隨機(jī)問(wèn)題，如何在試驗(yàn)階段對(duì)其進(jìn)行預(yù)測(cè)和提出解決方案呢？另一方面，某些系統(tǒng)在多次實(shí)驗(yàn)中要求相同的試驗(yàn)條件，如此才能更好地評(píng)估性能和發(fā)現(xiàn)問(wèn)題，這在實(shí)際環(huán)境下是不可能實(shí)現(xiàn)的，而系統(tǒng)仿真則可以對(duì)同一試驗(yàn)條件的精準(zhǔn)再現(xiàn)；其四，對(duì)于載人航天來(lái)說(shuō)，每一次發(fā)射都要耗費(fèi)漫長(zhǎng)的時(shí)間和高昂的金錢(qián)成本，次次都在真實(shí)的系統(tǒng)上進(jìn)行試驗(yàn)成本過(guò)高，系統(tǒng)仿真既節(jié)省了試驗(yàn)時(shí)間，又降低了試驗(yàn)成本。
仿真的基本方法是建立系統(tǒng)的結(jié)構(gòu)模型和量化分析模型，對(duì)模型進(jìn)行仿真試驗(yàn)。根據(jù)模型的不同，可以將系統(tǒng)仿真方法分為三類(lèi)：物理仿真、數(shù)學(xué)仿真和半實(shí)物仿真。物理仿真是按照真實(shí)系統(tǒng)的物理性質(zhì)構(gòu)造系統(tǒng)的物理模型，并在物理模型上進(jìn)行試驗(yàn)。物理仿真雖然直觀形象，但具有較多實(shí)驗(yàn)限制，且成本較高。數(shù)學(xué)仿真是對(duì)實(shí)際系統(tǒng)進(jìn)行抽象，并將其特性用數(shù)學(xué)關(guān)系加以描述而得到系統(tǒng)的數(shù)學(xué)模型，將其轉(zhuǎn)換為適合計(jì)算機(jī)編程的仿真模型，進(jìn)行仿真試驗(yàn)?，F(xiàn)代計(jì)算機(jī)技術(shù)的發(fā)展為系統(tǒng)仿真提供了條件。數(shù)學(xué)仿真在通信、控制、電力等多個(gè)工程和科研領(lǐng)域具有廣泛應(yīng)用。圖1、圖2是系統(tǒng)仿真的簡(jiǎn)單示例。
圖1 雙質(zhì)量-彈簧-阻尼系統(tǒng)
圖2 低通RC濾波器
由于連續(xù)系統(tǒng)和離散系統(tǒng)的數(shù)學(xué)模型具有很大差別，所以數(shù)學(xué)仿真可以分為連續(xù)系統(tǒng)仿真和離散系統(tǒng)仿真。半實(shí)物仿真則是對(duì)物理仿真和數(shù)學(xué)仿真的綜合應(yīng)用，即對(duì)系統(tǒng)中運(yùn)行規(guī)律明確的部分且簡(jiǎn)單的部分進(jìn)行數(shù)學(xué)仿真，規(guī)律不明確或難以數(shù)學(xué)建模的部分進(jìn)行物理建模，二者聯(lián)合起來(lái)完成仿真試驗(yàn)。例如在通信領(lǐng)域，對(duì)一個(gè)通信設(shè)備的系統(tǒng)仿真，除了通過(guò)數(shù)學(xué)建模完成對(duì)其通信系統(tǒng)的仿真分析，還要通過(guò)對(duì)設(shè)備使用環(huán)境進(jìn)行物理仿真，以確保設(shè)備在真實(shí)環(huán)境下的運(yùn)轉(zhuǎn)情況。
系統(tǒng)仿真是科研和工程領(lǐng)域不可缺少的環(huán)節(jié)，對(duì)工程技術(shù)和科學(xué)研究的發(fā)展具有重要意義。
[1] 福州大學(xué)《系統(tǒng)仿真技術(shù)》研究生課程
[2] http://www.bilibili.com/video/BV16E411A7zv?p=3&vd_source=f52c61024bca044c638d9c2fb92b3ede.
系統(tǒng)仿真軟件是一種用于模擬實(shí)際系統(tǒng)的軟件，它可以幫助我們更好地理解系統(tǒng)的運(yùn)行情況。仿真軟件是用于建模和分析復(fù)雜系統(tǒng)的強(qiáng)大工具。它被用于從工程到金融的各種行業(yè)，以幫助了解系統(tǒng)如何工作以及如何改進(jìn)。不同的仿真軟件有不同的特點(diǎn)，因此在選擇仿真軟件時(shí)，我們需要考慮不同的系統(tǒng)仿真軟件的優(yōu)缺點(diǎn)。
MATLAB/Simulink是一款功能強(qiáng)大的系統(tǒng)仿真軟件，它可以幫助我們快速構(gòu)建復(fù)雜的系統(tǒng)模型，并且可以模擬多種不同的系統(tǒng)。MATLAB是一個(gè)功能強(qiáng)大的通用工具，可用于模擬各種系統(tǒng)。它易于使用，具有廣泛的功能，包括圖形用戶(hù)界面、強(qiáng)大的編程語(yǔ)言和預(yù)構(gòu)建模型庫(kù)。MATLAB也是高度可擴(kuò)展的，允許用戶(hù)創(chuàng)建自己的模型和仿真。此外，MATLAB/Simulink還提供了豐富的可視化工具，可以幫助我們更好地理解系統(tǒng)的運(yùn)行情況。Simulink是一個(gè)圖形編程環(huán)境，允許用戶(hù)創(chuàng)建和模擬復(fù)雜系統(tǒng)。它被設(shè)計(jì)為易于使用，并具有廣泛的功能，包括一個(gè)預(yù)先構(gòu)建的模型庫(kù)、一種強(qiáng)大的編程語(yǔ)言和一個(gè)圖形用戶(hù)界面。Simulink還具有高度的可擴(kuò)展性，允許用戶(hù)創(chuàng)建自己的模型和仿真。但是，MATLAB/Simulink的上手較難，對(duì)于初學(xué)者來(lái)說(shuō)，可能需要花費(fèi)更多的時(shí)間來(lái)學(xué)習(xí)。
LabVIEW是另一款功能強(qiáng)大的系統(tǒng)仿真軟件，它可以幫助我們快速構(gòu)建復(fù)雜的系統(tǒng)模型，并且可以模擬多種不同的系統(tǒng)。此外，LabVIEW還提供了豐富的可視化工具，可以幫助我們更好地理解系統(tǒng)的運(yùn)行情況。LabVIEW主要用于開(kāi)發(fā)測(cè)量或控制系統(tǒng)。
Modelica是一種用于建模和仿真復(fù)雜系統(tǒng)的面向?qū)ο蟮木幊陶Z(yǔ)言。它提供了一種靈活的方法來(lái)描述系統(tǒng)的行為，并且可以用于模擬和分析復(fù)雜的系統(tǒng)。Modelica語(yǔ)言支持多種類(lèi)型的模型，包括動(dòng)態(tài)系統(tǒng)、靜態(tài)系統(tǒng)、熱系統(tǒng)、流體系統(tǒng)和電氣系統(tǒng)。它還支持多種類(lèi)型的仿真，包括時(shí)間域仿真、頻率域仿真和狀態(tài)空間仿真。
OpenModelica屬于開(kāi)源模擬軟件包。OpenModelica是一個(gè)功能強(qiáng)大的多功能工具，可用于模擬各種系統(tǒng)。它被設(shè)計(jì)為易于使用，并具有廣泛的功能，包括一個(gè)預(yù)先構(gòu)建的模型庫(kù)、一種強(qiáng)大的編程語(yǔ)言和一個(gè)圖形用戶(hù)界面。
Dymola基于Modelica語(yǔ)言提供一個(gè)基于模型的仿真環(huán)境，它可以用于模擬和分析復(fù)雜的系統(tǒng)。它提供了一個(gè)可視化的編程環(huán)境，可以幫助用戶(hù)快速構(gòu)建模型，并使用圖形化的工具進(jìn)行模擬和分析。與Matlab/Simulink不同，Dymola支持模型驅(qū)動(dòng)的仿真，可以更好地模擬復(fù)雜的系統(tǒng)。此外，Dymola還支持多種仿真技術(shù)，如動(dòng)態(tài)系統(tǒng)建模、模型驗(yàn)證和可視化等。
Scilab是一款開(kāi)源的科學(xué)計(jì)算軟件，它可以用于數(shù)值計(jì)算、矩陣運(yùn)算、繪圖、科學(xué)可視化、數(shù)據(jù)分析、模擬仿真等。它擁有豐富的函數(shù)庫(kù)，可以滿(mǎn)足各種科學(xué)計(jì)算需求。Scilab可以運(yùn)行在Windows、Linux、Mac OS X等操作系統(tǒng)上，支持多種編程語(yǔ)言，如C、C++、Fortran等，可以與其他軟件進(jìn)行交互，如Matlab、Maple、Mathematica等。
總體來(lái)說(shuō)，對(duì)于控制、通信、電氣等相關(guān)領(lǐng)域的理論研究推薦MATLAB/Simulink，LabVIEW則推薦用于測(cè)控系統(tǒng)，Modelica是一種編程語(yǔ)言，其自帶的求解器較為弱，偏好該語(yǔ)言的用戶(hù)可以嘗試Dymola等相關(guān)軟件，Scilab則是一個(gè)較為小眾的開(kāi)源軟件。
一、 最小二乘原理介紹
最小二乘法是一種常見(jiàn)的數(shù)學(xué)方法，用于解決線(xiàn)性回歸問(wèn)題。它可以用來(lái)估計(jì)因變量（目標(biāo)變量）與自變量之間的線(xiàn)性關(guān)系，即通過(guò)已知的一組數(shù)據(jù)來(lái)確定未知的回歸方程。
具體來(lái)說(shuō)，最小二乘法的目標(biāo)是找到一條直線(xiàn)（或者更一般的曲線(xiàn)），使得這條直線(xiàn)與已知數(shù)據(jù)點(diǎn)的距離平方和最小。這條直線(xiàn)的方程稱(chēng)為回歸方程，可以用來(lái)預(yù)測(cè)因變量的值。
最小二乘法的原理是通過(guò)最小化誤差平方和來(lái)確定回歸方程的系數(shù)。誤差指的是已知數(shù)據(jù)點(diǎn)到回歸直線(xiàn)的距離，也就是因變量的真實(shí)值與預(yù)測(cè)值之間的差異。通過(guò)最小化誤差平方和，可以找到最優(yōu)的回歸系數(shù)，使得預(yù)測(cè)值與真實(shí)值之間的誤差最小。
在實(shí)際應(yīng)用中，最小二乘法常常用于數(shù)據(jù)擬合和預(yù)測(cè)。例如，在金融領(lǐng)域，最小二乘法可以用來(lái)預(yù)測(cè)股票價(jià)格和交易量；在工程領(lǐng)域，最小二乘法可以用來(lái)擬合實(shí)驗(yàn)數(shù)據(jù)和優(yōu)化設(shè)計(jì)參數(shù)。
在計(jì)算中，最小二乘法可以通過(guò)矩陣運(yùn)算來(lái)求解。具體來(lái)說(shuō)，可以使用線(xiàn)性代數(shù)中的矩陣求解方法，例如求解矩陣的逆或者使用QR分解。在北太天元中，可以使用polyfit函數(shù)來(lái)實(shí)現(xiàn)最小二乘法求解線(xiàn)性回歸問(wèn)題。
二、 根據(jù)最小二乘原理進(jìn)行多項(xiàng)式擬合
原始數(shù)據(jù)點(diǎn)為，根據(jù)最小二乘法進(jìn)行多項(xiàng)式擬合：
(1) 設(shè)擬合多項(xiàng)式為
(2) 各點(diǎn)到這條曲線(xiàn)的距離之和，即誤差平方和為：
(3) 根據(jù)最小二乘原理，等式右邊求ai偏導(dǎo)數(shù)：
(4) 化簡(jiǎn)得：
(5) 化成矩陣形式：
(6) 化簡(jiǎn)得
(7) 求解上述方程組，即可得到系數(shù)矩陣，從而得到擬合多項(xiàng)式。
三、 基于北太天元實(shí)現(xiàn)最小二乘法，并對(duì)比擬合結(jié)果與polyfit結(jié)果
附錄代碼：
clear all
close all
clc
N=10;%多項(xiàng)式階數(shù)
%原始數(shù)據(jù)點(diǎn)
X=10:0.05:15;
M=length(X);
p=[5.6,2.5,3.3,0.18,0.29,2.4,4.8,1.2,7.1,1];
Y=zeros(1,M);
for i=1:10
Y=Y+p(i)*X.^(i-1);
end
% Y=Y+unifrnd(-2e9,2e9,1,M);
%多項(xiàng)式擬合
X1=zeros(N+1,M);
for i=1:M
X1(1,i)=1;
end
for i=2:N+1
X1(i,:)=power(X,i-1); %構(gòu)造系數(shù)矩陣
end
XX=X1';
P_n=(XX'*XX)\XX'*Y';
X_n=min(X):max(X);
Y_n=zeros(1,length(X_n));
for i=1:N+1
Y_n=Y_n+P_n(i)*X_n.^(i-1);
end
%polyfit
P=zeros(1,N);
P=polyfit(X,Y,N);
X_polyfit=zeros(1,M);
Y_polyfit=zeros(1,M);
X_polyfit=linspace(min(X),max(X),M);
Y_polyfit=polyval(P,X_polyfit);
figure(1);
plot(X,Y,'.');
hold on
plot(X_n,Y_n,'g');
hold on
plot(X_polyfit,Y_polyfit,'r');
legend('原始數(shù)據(jù)點(diǎn)','多項(xiàng)式擬合','polyfit','Location','northwest');
該求解器有變步長(zhǎng)和定步長(zhǎng)兩種類(lèi)型。不同類(lèi)型有著不同的求解器，其中ode45求解器屬于變步長(zhǎng)的一種，采用Runge-Kutta算法。
ode45表示采用四階-五階Runge-Kutta算法，它用4階方法提供候選解，5階方法控制誤差，是一種自適應(yīng)步長(zhǎng)（變步長(zhǎng)）的常微分方程數(shù)值解法，其整體截?cái)嗾`差為(Δx)^5。解決的是非剛性常微分方程。
ode45是解決數(shù)值解問(wèn)題的首選方法，若長(zhǎng)時(shí)間沒(méi)結(jié)果，應(yīng)該就是剛性的，可換用ode15s試試。
在各種龍格－庫(kù)塔法當(dāng)中有一個(gè)方法十分常用，以至于經(jīng)常被稱(chēng)為“RK4”或者就是“龍格－庫(kù)塔法”。該方法主要是在已知方程導(dǎo)數(shù)和初值信息，利用計(jì)算機(jī)仿真時(shí)應(yīng)用，省去求解微分方程的復(fù)雜過(guò)程。
令初值問(wèn)題表述如下。
y' =f(t,y)，y(t0)=y0
則，對(duì)于該問(wèn)題的RK4由如下方程給出：
y(n+1)=y(n)+h/6*(k1+2*k2+2*k3+k4)
其中
k1=f(t(n),y(n))
k2=f(t(n)+h/2,y(n)+h/2*k1)
k3=f(t(n)+h/2,y(n)+h/2*k2)
k4=f(t(n)+h,y(n)+h*k3)
這樣，下一個(gè)值（yn+1）由現(xiàn)在的值（yn）加上時(shí)間間隔（h）和一個(gè)估算的斜率的乘積所決定。該斜率是以下斜率的加權(quán)平均：
k1是時(shí)間段開(kāi)始時(shí)的斜率；
k2是時(shí)間段中點(diǎn)的斜率，通過(guò)歐拉法采用斜率k1來(lái)決定y在點(diǎn)tn+h/2的值；
k3也是中點(diǎn)的斜率，但是這次采用斜率k2決定y值；
k4是時(shí)間段終點(diǎn)的斜率，其y值用k3決定。
當(dāng)四個(gè)斜率取平均時(shí)，中點(diǎn)的斜率有更大的權(quán)值：
RK4法是四階方法，也就是說(shuō)每步的誤差是h階，而總積累誤差為h階。
求解微分方程d2y/dt2- 7(1-y2)dy/dt + y = 0在初始條件下y(0) = 1、y’(0) = 0下的解，并畫(huà)出解的圖。
解：設(shè)x1 = y，x2 = dy/dt，將二階方程化成一階方程組
dx1 = x2 x1(0) = 1;
dx1/dt = 7(1-x12)*x2-x1 x2(0) = 0;
在北太天元腳本編輯器窗口輸入以下程序，并保存文件名為vdp.m。
function fy = vdp(t,x)
fy = [x(2);7*(1-x(1)^2)*x(2)-x(1)];
在北太天元命令行窗口輸入以下程序。
>> y0 = [1;0];
>> [t,x] = ode45(@vdp,[0 40],y0);
>> y = x(:,1);
>> plot(t,y)
>> xlabel('t')
>> ylabel('y')
輸出結(jié)果如下圖所示：
用雷達(dá)跟蹤目標(biāo)，目標(biāo)的運(yùn)動(dòng)可以看成是在徑向和橫向內(nèi)的二維運(yùn)動(dòng)，其運(yùn)動(dòng)方程和觀測(cè)方程分別為：
s1(k)、v1(k)和y1(k)分別為徑向距離、速度和觀測(cè)值，而s2(k)、v2(k)和y2(k)分別為橫向距離、速度和觀測(cè)值。u1(k)和u2(k)是狀態(tài)噪聲，是目標(biāo)速度的波動(dòng)；w1(k)和w2(k)是觀測(cè)噪聲；四種噪聲的均值都為0，呈高斯分布，互不相關(guān)。
狀態(tài)方程和輸出方程：
其中：
在推導(dǎo)上式時(shí)，wi(k)和wi(k+1)是隨機(jī)過(guò)程中不同時(shí)刻的兩個(gè)隨機(jī)變量，我們認(rèn)為這兩個(gè)隨機(jī)變量統(tǒng)計(jì)獨(dú)立，而且w(k)是平穩(wěn)隨機(jī)過(guò)程，其不同時(shí)刻的方差相同。兩個(gè)時(shí)刻的時(shí)差T為雷達(dá)掃描一圈的時(shí)間。
仿真結(jié)果如下：
clear allclose allclcT=1;sigma_u1=300;sigma_u2=0.12;sigma_w1=900000;sigma_w2=900000;N=500;u1=sqrt(sigma_u1).*randn(1,N);u2=sqrt(sigma_u2).*randn(1,N);w1=sqrt(sigma_w1).*randn(1,N);w2=sqrt(sigma_w2).*randn(1,N);A=[ 1,T,0,0; 0,1,0,0; 0,0,1,T; 0,0,0,1; ];C=[ 1,0,0,0; 0,0,1,0; ];Q=[ %狀態(tài)噪聲協(xié)方差矩陣 0,0,0,0; 0,sigma_u1,0,0; 0,0,0,0; 0,0,0,sigma_u2; ];R=[ %觀測(cè)噪聲協(xié)方差矩陣 sigma_w1,0; 0,sigma_w2; ];I=[ 1,0,0,0; 0,1,0,0; 0,0,1,0; 0,0,0,1; ];x0=[10000;300;900;256;];%初始值% y0=[10000;0;];x=zeros(4,N);y=zeros(2,N);x(:,1)=x0;% y(:,1)=y0;x_calculation=zeros(4,N); %狀態(tài)濾波估計(jì)值u=[zeros(1,N);u1;zeros(1,N);u2];w=[w1;w2];y(:,1)=C*x(:,1)+w(:,1);P2=[ sigma_w1, sigma_w1/T, 0, 0; sigma_w1/T,sigma_u1+2*sigma_w1/(T^2),0, 0; 0, 0, sigma_w2, sigma_w2/T; 0, 0, sigma_w2/T,sigma_u2+2*sigma_w2/(T^2); ];p=P2;for i=1:N-1 x(:,i+1)=A*x(:,i)+u(:,i);%真實(shí)值 y(:,i+1)=C*x(:,i+1)+w(:,i+1); endx_calculation(:,2)=[y(1,1) , (y(1,2)-y(1,1))/T ,y(2,2), (y(2,2)-y(2,1))/T]';for i=3:N x1=A*x_calculation(:,i-1);%狀態(tài)預(yù)測(cè)方程 p1=A*p*A'+Q; %狀態(tài)預(yù)測(cè)誤差協(xié)方差矩陣 K=p1*C'*inv(C*p1*C'+R); %最佳增益方程 x_calculation(:,i)=x1+K*(y(:,i)-C*x1); %濾波估值方程 p=(I-K*C)*p1; %濾波估值誤差 協(xié)方差方程end t=1:N; figure(4)plot(t,y(1,:),'b',t,x_calculation(1,:),'r');xlabel('t');ylabel('s1');legend('徑向距離觀測(cè)值','徑向距離估計(jì)值','Location','southeast');title('二維卡爾曼濾波器');figure(5)plot(t,y(2,:),'b',t,x_calculation(3,:),'r');xlabel('t');ylabel('s2');legend('橫向距離觀測(cè)值','橫向距離估計(jì)值','Location','southeast');1) 忽略轉(zhuǎn)向系的影響，以前、后輪轉(zhuǎn)角作為輸入；
2) 汽車(chē)只進(jìn)行平行于地面的平面運(yùn)動(dòng)，而忽略懸架的作用；
3) 汽車(chē)前進(jìn)（縱軸）速度不變，只有沿y軸的側(cè)向速度和繞z軸的橫擺運(yùn)動(dòng)(ay<0.4g) ；
4) 驅(qū)動(dòng)力不大，對(duì)側(cè)偏特性無(wú)影響；
5) 忽略空氣阻力；
6) 忽略左右輪胎因載荷變化引起輪胎特性的變化；
7) 忽略回正力矩的變化。
根據(jù)模型假設(shè)建立如圖1所示的二自由度汽車(chē)模型。
對(duì)模型受力分析，存在3個(gè)方向的受力平衡，分別為x、y和繞Z的力矩平衡，建立力學(xué)方程如下。
在baltamulink中搭建狀態(tài)空間模型，模型如圖所示。
（1）在前輪偏轉(zhuǎn)角為1，后輪偏轉(zhuǎn)角為1，車(chē)速為40km/h的情況下，輸出前后輪的橫向位移情況，輸出結(jié)果如圖3。
圖3
通過(guò)建立汽車(chē)動(dòng)力學(xué)模型，對(duì)汽車(chē)操縱性進(jìn)行餓模擬。根據(jù)仿真結(jié)果可以發(fā)現(xiàn)車(chē)速和前輪轉(zhuǎn)角都對(duì)二自由度汽車(chē)的操縱穩(wěn)定性有很大影響。汽車(chē)以較低速度、較小的前輪轉(zhuǎn)角行駛時(shí)，是相對(duì)安全的。
通過(guò)分析圖3可以看出前、后輪的橫向位移都是發(fā)散的，這是因?yàn)榻o前輪的一個(gè)階躍響應(yīng)，一直存在前輪轉(zhuǎn)角，同時(shí)系統(tǒng)沒(méi)有加入閉環(huán)控制，屬于開(kāi)環(huán)控制，這就導(dǎo)致前后輪的橫向位移處于發(fā)散狀態(tài)。
baltamatica代碼如下
clc;clear;close all;
%% 基本車(chē)輛參數(shù)
v=40/3.6;%輸入為km/h，方程單位為m/s
m=16000;%車(chē)重
I=10.85*m;%轉(zhuǎn)動(dòng)慣量
cf=340000;%側(cè)偏剛度
cr=cf;
lf=2.65;%前軸到重心的距離
lr=3.35;%后軸到重心的距離
a=pi/180;
%% 組裝矩陣
A=[-(cf+cr)/(m*v) -1-(cf*lf-cr*lr)/(m*v^2) 0 0
-(cf*lf-cr*lr)/I -(cf*lf^2+cr*lr^2)/(I*v) 0 0
v lf 0 0
v -lr 0 0];
B=[cf/(m*v) cr/(m*v)
cf*lf/I cr*lr/I
0 0
0 0];
C=[0 0 1 0
0 0 0 1];
D=zeros(2,2);
問(wèn)題：?jiǎn)钨|(zhì)量彈簧阻尼機(jī)械系統(tǒng)類(lèi)似于下圖，外力f(x)為輸入量；質(zhì)量位移x(t)為輸出量；質(zhì)量m為1kg；剛度為5N/m；阻尼系數(shù)f為0.3N·s/m。繪制系統(tǒng)位移輸出響應(yīng)曲線(xiàn)。
圖1：?jiǎn)钨|(zhì)量彈簧阻尼機(jī)械系統(tǒng)
首先：?jiǎn)钨|(zhì)量彈簧阻尼機(jī)械系統(tǒng)的震動(dòng)方程為：
m dx2/dt2 + c dx/dt + kx = f(x)
取狀態(tài)向量為X(t) = [x(t) dx/dt]’，輸出向量為U=[f(x)]，輸出向量為Y=[x(t)]，則狀態(tài)方程為：
dX/dt = AX(t) + BU
Y = CX(t) + DU
式中，A = [0 1; -k/m -c/m]; B = [0; 1/m]; C = [1 0]; D = 0;(m = 1; k = 5; c = 0.3)。
其次，通過(guò)計(jì)算得到傳遞函數(shù)的分子分母系數(shù)為：[1]； [1 0.3 5];
即，傳遞函數(shù)為
G(s) = 1/(s^2 + 0.3 s + 5);
最后，利用北太真元建立傳遞函數(shù)仿真模型，如下圖所示：
設(shè)置：仿真時(shí)長(zhǎng)：10s
求解器為：定步長(zhǎng) ode4；
步長(zhǎng)：1s
輸入常量10
仿真結(jié)果如下圖所示：
在MATLAB中創(chuàng)建相同模型，仿真時(shí)間、仿真求解器、步長(zhǎng)均相同；仿真模型和仿真結(jié)果如下圖所示：
通過(guò)對(duì)比發(fā)現(xiàn)，北太真元計(jì)算結(jié)果與MATLAB仿真結(jié)果完全一致。
問(wèn)題：
二階電路動(dòng)態(tài)系統(tǒng)中，該系統(tǒng)是由電阻R、電感L和電容C組成的無(wú)源網(wǎng)絡(luò)；ui(t)為輸入，i(t)為電流，u0(t)為輸出；設(shè)R=1Ω，L= 2H，C=2F；系統(tǒng)的初始狀態(tài)為0，外加的輸入為單位階躍信號(hào)。求系統(tǒng)的輸出波形。
圖：二階電路動(dòng)態(tài)系統(tǒng)
首先，在解決該問(wèn)題時(shí)，需要了解基爾霍夫電壓定律；它的內(nèi)容是，在任何一個(gè)閉合回路中，各元件上的電壓降的代數(shù)和等于電動(dòng)勢(shì)的代數(shù)和，即從一點(diǎn)出發(fā)繞回路一周回到該點(diǎn)時(shí)，各段電壓的代數(shù)和恒等于零，即∑U=0。
基爾霍夫電壓定律表明：
沿著閉合回路所有元件兩端的電勢(shì)差（電壓）的代數(shù)和等于零。
或者描述為：
沿著閉合回路的所有電動(dòng)勢(shì)的代數(shù)和等于所有電壓降的代數(shù)和。
以方程表達(dá)，對(duì)于電路的任意閉合回路有：
其中，m 是這閉合回路的元件數(shù)目，vk 是元件兩端的電壓，可以是實(shí)數(shù)或復(fù)數(shù)。
基爾霍夫電壓定律不僅應(yīng)用于閉合回路，也可以把它推廣應(yīng)用于回路的部分電路。
根據(jù)上述原理，可寫(xiě)出上面問(wèn)題的回路方程如下：
L di(t)/dt + 1/C∫i(t)dt + Ri(t) = ui(t)
消去中間變量i(t)，可以得到描述網(wǎng)絡(luò)輸入輸出關(guān)系的微分方程為：
LC d2u0(t)/dt2 + RC du0(t)/dt + u0(t) = ui(t)
帶入電路參數(shù)R=1Ω，L= 2H，C=2F；整理可得：
d2u0(t)/dt2 + 0.5 du0(t)/dt + 0.25u0(t) = 0.25ui(t)
從方面方程可以得到，傳遞函數(shù)分子等式右邊系數(shù)：[0.25]；分母為左邊等式系數(shù)：[1 0.5 0.25]
根據(jù)該系數(shù)，在北太真元設(shè)計(jì)傳遞函數(shù)仿真模型，再加入一個(gè)階躍信號(hào)模塊，得到仿真模型如下圖所示：
設(shè)置：仿真時(shí)長(zhǎng)：40s
求解器為：定步長(zhǎng)；
步長(zhǎng)：1s；
得到仿真結(jié)果如下圖所示：
在MATLAB中創(chuàng)建相同模型，仿真時(shí)間、仿真求解器、步長(zhǎng)均相同；仿真模型和仿真結(jié)果如下圖所示：
通過(guò)對(duì)比發(fā)現(xiàn)，北太真元計(jì)算結(jié)果與MATLAB仿真結(jié)果完全一致。
由于汽車(chē)在行走時(shí)，路面不平，汽車(chē)行駛中的路面可簡(jiǎn)化成正弦函數(shù)。
可把汽車(chē)行走的路面看做激勵(lì)，忽略輪胎的彈性與質(zhì)量，得到分析車(chē)身垂直振動(dòng)的最簡(jiǎn)單的單質(zhì)量系統(tǒng)，適用于低頻激勵(lì)情況。汽車(chē)行駛可看作如下模型：
上圖為單一自由度系統(tǒng)的簡(jiǎn)圖，設(shè)X(t)及Xs(t)分別是質(zhì)量塊及支承的位移，支承的運(yùn)動(dòng)規(guī)律是：
Xs =asinwt
由于支承的運(yùn)動(dòng)，質(zhì)量塊收到的彈性恢復(fù)力為k(X - Xs)，阻尼力為c(Vx-Vxs)
根據(jù)達(dá)朗伯原理可得如下的運(yùn)動(dòng)微分方程：
由（1）和（2）得：
在此系統(tǒng)中除了有彈性恢復(fù)力及阻尼力作用外，還始終作用于簡(jiǎn)諧激勵(lì)力：
Px=P0sinwt
簡(jiǎn)諧激勵(lì)：
激勵(lì)隨時(shí)間的變化規(guī)律可用正弦或余弦函數(shù)表示；
振動(dòng)響應(yīng)亦為時(shí)間的正弦和余弦函數(shù)（簡(jiǎn)諧振動(dòng)）。
結(jié)合上面的運(yùn)動(dòng)微分方程和簡(jiǎn)諧激勵(lì)力方程，可得系統(tǒng)的運(yùn)動(dòng)微分方程為：
令：物體質(zhì)：m = 1 kg，彈簧剛度：k = 3 N/m，阻尼：c = 4 N·s/m，作用力P = 2sin(2t + π/3)，研究物體的位移隨時(shí)間的變化規(guī)律。
通過(guò)北太真元建立系統(tǒng)運(yùn)動(dòng)微分方程模型，如下圖所示：
設(shè)置參數(shù)：
正弦波產(chǎn)生模塊：幅值：2；偏置：0；頻率(弧度/秒)：2；相位(弧度)：pi/3≈1；
一階積分模塊：積分初始值：0.5；
增益模塊：增益數(shù)值：4；
常量模塊：常量值：3；
結(jié)束時(shí)間：10s；
求解器：ode4；
步長(zhǎng)：0.01s。
得到的仿真結(jié)果，如下圖所示：
問(wèn)題：利用汽車(chē)橫擺角速度傳遞函數(shù)和質(zhì)心側(cè)偏角傳遞函數(shù)，對(duì)汽車(chē)時(shí)域響應(yīng)進(jìn)行仿真，繪制汽車(chē)橫擺角速度和質(zhì)心偏側(cè)角的時(shí)域特性曲線(xiàn)。汽車(chē)時(shí)域響應(yīng)仿真所需參數(shù)見(jiàn)下表。
由于汽車(chē)橫擺角速度傳遞函數(shù)（G1(s)）和質(zhì)心偏側(cè)角傳遞函數(shù)（G2(s)）分別為：
G1(s) = ((s - a11)*b21 + a21*b11) / s2- (a11 + a22)*s + a11*a22 - a12*a21;
G2(s) = ((s - a22)*b11 + a12*b21) / s2- (a11 + a22)*s + a11*a22 - a12*a21;
式中，
a11 = (Ka1 + Ka2) / mu; a12 = (aKa1 - bKa2 - mu2) / mu2;
a21 = (aKa1 - bKa2) / Iz; a22 = (a2Ka1 + b2Ka2) / Iz*u;
b11 = -Ka1 / mu; b21 = -aKa1/Iz；
汽車(chē)速度分別選取10m/s、20m/s、30m/s；在仿真時(shí)間0s時(shí)給前輪一個(gè)階躍信號(hào)，使前輪轉(zhuǎn)角從0°轉(zhuǎn)到15°，并保持不變。根據(jù)汽車(chē)橫擺角速度傳遞函數(shù)和質(zhì)心偏移角傳遞函數(shù)，建立模型，繪制不同車(chē)速下的汽車(chē)橫擺角度和質(zhì)心側(cè)偏角的時(shí)域特性曲線(xiàn)。
首先：通過(guò)北太天元計(jì)算汽車(chē)橫擺角度速度傳遞函數(shù)分子和分母的系數(shù)，在北太天元依次輸入下面語(yǔ)句；
>> m=3018; Iz=10437; a=1.84; b=1.88; k1=-23147; k2=-38318;
>> u= [10 20 30];
>> a11 = (k1 + k2)/m./u;
>> a12 = (a*k1 - b*k2 -m.*u^2)/m./u^2;
>> a21 = (a * k1 - b * k2)/Iz;
>> a22=(a^2*k1 + b^2*k2)/Iz./u;
>> b11 = -k1/m./u;
>> b21 = -a*k1/Iz;
>> b1 = b21;
>> b2 = a21*b11-a11*b21;
>> b3 = -a11-a22;
>> b4 = a11.*a22-a12.*a21;
>> num = [b1 b2];
>> den = [1,b3,b4];
得到結(jié)果如下圖1所示 ；將命令行窗口，和工作區(qū)窗口放大后如圖2、圖3所示；
圖1
圖2
圖3
因?yàn)?，汽?chē)橫擺角度速度傳遞函數(shù)的分子系數(shù)為：num = [b1, b2]; 分母系數(shù)為：den = [1, b3, b4]; 所以，從圖3紅色框中可以得到各項(xiàng)系數(shù)如下：
當(dāng)汽車(chē)速度 s = 10 m/s 時(shí)，分子系數(shù)：num = [4.0807 10.4748]; 分母系數(shù)：den = [1 4.0851 6.7181];
當(dāng)汽車(chē)速度 s = 20 m/s 時(shí)， 分子系數(shù)：num = [4.0807 5.2347]; 分母系數(shù)：den = [1 2.0425 3.7956];
當(dāng)汽車(chē)速度 s = 30 m/s 時(shí)， 分子系數(shù)：num = [4.0807 3.4916]; 分母系數(shù)：den = [1 1.3617 3.2544];
又因?yàn)椋诜抡鏁r(shí)間0s時(shí)給前輪一個(gè)階躍信號(hào)，使前輪轉(zhuǎn)角從0°轉(zhuǎn)到15°；所以模型還需一個(gè)階躍信號(hào)模塊，階躍時(shí)間=0；且，還需一個(gè)增益模塊，增益= pi*15/180 = 0.2618。
通過(guò)北太真元建立汽車(chē)橫擺角度速度傳遞函數(shù)模型，如下圖所示：
設(shè)置參數(shù)：
仿真時(shí)長(zhǎng)：10s；步長(zhǎng)0.01s；求解器：ode4
得到的仿真結(jié)果，如下圖所示：
紫色代表速度10m/s時(shí)響應(yīng)曲線(xiàn)；
墨綠色代表速度20m/s時(shí)響應(yīng)曲線(xiàn)；
橙色代表速度30m/s時(shí)響應(yīng)曲線(xiàn)；
首先：通過(guò)北太天元計(jì)算質(zhì)心偏側(cè)角傳遞函數(shù)分子和分母的系數(shù)，在北太天元依次輸入下面語(yǔ)句；
>> m=3018;Iz=10437;a=1.84;b=1.88;k1=-23147;k2=-38318;
>> u= [10 20 30];
>> a11 = (k1 + k2)/m./u;
>> a12 = (a*k1 - b*k2 -m.*u.^2)/m./u.^2;
>> a21 = (a * k1 - b * k2)/Iz;
>> a22=(a^2*k1 + b^2*k2)/Iz./u;
>> b11 = -k1/m./u;
>> b21 = -a*k1/Iz;
>> b1 = b11;
>> b2 = a12*b21-a22*b11;
>> b3 = -a11-a22;
>> b4 = a11.*a22-a12.*a21;
>> num = [b1 b2];
>> den = [1,b3,b4];
命令行窗口如圖4所示；參數(shù)計(jì)算結(jié)果如圖5所示；
圖4
圖5
因?yàn)?，質(zhì)心偏側(cè)角傳遞函數(shù)的分子系數(shù)為：num = [b1, b2]; 分母系數(shù)為：den = [1, b3, b4]; 所以，從圖5紅色框中可以得到各項(xiàng)系數(shù)如下：
當(dāng)汽車(chē)速度 s = 10 m/s 時(shí)，分子系數(shù)：num = [0.766965 -2.11146]; 分母系數(shù)：den = [1 4.08507 6.71806];
當(dāng)汽車(chē)速度 s = 20 m/s 時(shí)， 分子系數(shù)：num = [0.383482 -3.58841]; 分母系數(shù)：den = [1 2.04254 3.7956];
當(dāng)汽車(chē)速度 s = 30 m/s 時(shí)， 分子系數(shù)：num = [0.255655 -3.86191]; 分母系數(shù)：den = [1 1.36169 3.2544];
又因?yàn)椋诜抡鏁r(shí)間0s時(shí)給前輪一個(gè)階躍信號(hào)，使前輪轉(zhuǎn)角從0°轉(zhuǎn)到15°；所以模型還需一個(gè)階躍信號(hào)模塊，階躍時(shí)間=0；且，還需一個(gè)增益模塊，增益= pi*15/180 = 0.2618。
通過(guò)北太真元建立質(zhì)心偏側(cè)角傳遞函數(shù)模型 同 汽車(chē)橫擺角度速度傳遞函數(shù)模型。
設(shè)置參數(shù)：
仿真時(shí)長(zhǎng)：10s；步長(zhǎng)0.01s；求解器：ode4
得到的仿真結(jié)果，如下圖所示：
淺綠色表速度10m/s時(shí)響應(yīng)曲線(xiàn)；
紫代表速度20m/s時(shí)響應(yīng)曲線(xiàn)；
墨綠色代表速度30m/s時(shí)響應(yīng)曲線(xiàn)；
實(shí)例：已知傳遞函數(shù)G(s) = ω2/ s2+ 2ζωs + ω2，分析阻尼系數(shù)和固有頻率對(duì)性能的影響。
（1）假設(shè) ω = 1, ζ = 0, 0.8, 1.5;
（2）假設(shè) ζ = 1, ω = 1, 2, 3;
從（1）可得，阻尼系數(shù)傳遞函數(shù)的系數(shù)可以是：
G(s1) = [1; 1 0 1];
G(s2) = [1; 1 1.6 1];
G(s3) = [1; 1 3 1];
通過(guò)北太真元建立“阻尼系數(shù)對(duì)系統(tǒng)性能的影響”模型，如下圖所示：
設(shè)置參數(shù)：仿真時(shí)長(zhǎng)：20s；步長(zhǎng)0.1s；求解器：ode4
得到的仿真結(jié)果，如下圖所示：
墨綠色代表阻尼系數(shù)=0時(shí)的特性曲線(xiàn)；
淡綠色代表阻尼系數(shù)=0.8時(shí)的特性曲線(xiàn)；
紅色代表阻尼系數(shù)=1.5時(shí)的特性曲線(xiàn)。
從結(jié)果圖中可以看出，阻尼系數(shù)決定了系統(tǒng)的振蕩幅度，阻尼系數(shù)越小，振蕩幅度越大。
從（2）可得，固有頻率傳遞函數(shù)的系數(shù)可以是：
G(s1) = [1; 1 0.5 1];
G(s2) = [4; 1 1 4];
G(s3) = [9; 1 1.5 9];
通過(guò)北太真元建立“固有頻率對(duì)系統(tǒng)性能的影響”模型，同上模型；
設(shè)置參數(shù)：仿真時(shí)長(zhǎng)：20s；步長(zhǎng)0.1s；求解器：ode4
得到的仿真結(jié)果，如下圖所示：
墨綠色代表固有頻率=1時(shí)的特性曲線(xiàn)；
淡綠色代表固有頻率=2時(shí)的特性曲線(xiàn)；
紫色代表固有頻率=3時(shí)的特性曲線(xiàn)。
從結(jié)果圖中可以看出，固有頻率決定了系統(tǒng)的振蕩頻率，固有頻率越大，系統(tǒng)的振蕩越高，響應(yīng)速度也越快。
問(wèn)題：利用汽車(chē)橫擺角速度傳遞函數(shù)和質(zhì)心側(cè)偏角傳遞函數(shù)，對(duì)汽車(chē)時(shí)域響應(yīng)進(jìn)行仿真，繪制汽車(chē)橫擺角速度和質(zhì)心偏側(cè)角的時(shí)域特性曲線(xiàn)。汽車(chē)時(shí)域響應(yīng)仿真所需參數(shù)見(jiàn)下表。
取狀態(tài)向量為X = [β ωr]’，輸入向量U = [δ1]，輸出向量為Y = [β ωr]’，狀態(tài)空間方程為：
式中，A = [(K1+K2)/mu, (aK1 - BK2)/mu2-1; (aK1-bK2)/Iz, (a2K1+b2K2)/Iz*u] 稱(chēng)為系統(tǒng)矩陣；B = [-K1/mu; -aK1/Iz] 稱(chēng)為控制矩陣；C = [1 0; 0 1] 稱(chēng)為輸出矩陣；D = [0; 0] 稱(chēng)為傳遞矩陣。
汽車(chē)速度分別選取20m/s、30m/s、40m/s；在仿真時(shí)間0s時(shí)給前輪一個(gè)階躍信號(hào)，使前輪轉(zhuǎn)角從0°轉(zhuǎn)到10°，并保持不變。根據(jù)汽車(chē)狀態(tài)空間模型，建立模型，繪制不同車(chē)速下的汽車(chē)橫擺角速度和質(zhì)心側(cè)偏角的時(shí)域特性曲線(xiàn)。
首先：通過(guò)北太天元計(jì)算汽車(chē)狀態(tài)空間方程的系統(tǒng)矩陣和控制矩陣，在北太天元依次輸入下面語(yǔ)句；
>> m=2050;Iz=5600;a=1.5;b=1.8;L=3.3;
>> k1=-38900;k2=-39200;
>> u= [20 30 40];
>> a11 = (k1 + k2)/m./u;
>> a12 = (a*k1 - b*k2 -m.*u.^2)/m./u.^2;
>> a21 = (a*k1 - b*k2)/Iz;
>> a22=(a^2*k1 + b^2*k2)/Iz./u;
>> b11 = -k1/m./u;
>> b21 = -a*k1/Iz;
得到結(jié)果如下圖1所示 ；
圖1
將命令行窗口，和工作區(qū)窗口放大后如圖2、圖3所示；
圖2
圖3
因?yàn)?，汽?chē)狀態(tài)空間方程的系統(tǒng)矩陣為：A = [a11, a12; a21, a22]；控制矩陣為：B = [ b11; b21]； 所以，從圖3紅色框中可以得到各項(xiàng)系數(shù)如下：
當(dāng)汽車(chē)速度 s = 20 m/s 時(shí)，系統(tǒng)矩陣：A = [-1.90488,-0.98511;2.18036,-1.91547]; 控制矩陣：B = [0.94878; 10.4196];
當(dāng)汽車(chē)速度 s = 30 m/s 時(shí)，系統(tǒng)矩陣：A = [-1.26992,-0.993382;2.180436,-1.27698]; 控制矩陣：B = [0.63252; 10.4196];
當(dāng)汽車(chē)速度 s = 40 m/s 時(shí)，系統(tǒng)矩陣：A = [-0.952439,-0.996277;2.180436,-0.957737]; 控制矩陣：B = [0.47439; 10.4196];
狀態(tài)方程輸出矩陣C = [1 0; 0 1]；傳遞矩陣D = [0; 0]。
又因?yàn)?，在仿真時(shí)間0s時(shí)給前輪一個(gè)階躍信號(hào)，使前輪轉(zhuǎn)角從0°轉(zhuǎn)到10°；所以模型還需一個(gè)階躍信號(hào)模塊，階躍時(shí)間=0；且，還需一個(gè)增益模塊，增益= pi*10/180 = 0.1745。
通過(guò)北太真元建立汽車(chē)狀態(tài)空間模型，如下圖所示：
設(shè)置參數(shù)：仿真時(shí)長(zhǎng)：10s；步長(zhǎng)0.01s；求解器：ode4
得到的仿真結(jié)果，如下圖所示：
上半部分代表汽車(chē)橫擺角速度時(shí)域特性曲線(xiàn)；即：
墨綠色代表速度20m/s時(shí)的特性曲線(xiàn)；
綠色代表速度30m/s時(shí)的特性曲線(xiàn)；
紅色代表速度40m/s時(shí)的特性曲線(xiàn)。
下半部分代表汽車(chē)質(zhì)心側(cè)偏角時(shí)域特性曲線(xiàn)；即：
紫色代表速度20m/s時(shí)的特性曲線(xiàn)；
橙色代表速度30m/s時(shí)的特性曲線(xiàn)；
藍(lán)色代表速度40m/s時(shí)的特性曲線(xiàn)。
使用4:1衰減曲線(xiàn)法設(shè)計(jì)下列被控傳遞函數(shù)的PI控制器，分別計(jì)算P控制、PI控制的參數(shù)值，并繪制控制前后系統(tǒng)的單位階躍響應(yīng)曲線(xiàn)。4:1衰減法控制參數(shù)計(jì)算公式如下表所示：
4:1衰減法控制
被控傳遞函數(shù)方程如下：
Gp(s) =1 / 100^3 + 80^s + 17s + -1;
調(diào)節(jié)參數(shù)時(shí)，比例系數(shù)由小變大，并增加擾動(dòng)觀察響應(yīng)過(guò)程，知道響應(yīng)曲線(xiàn)峰值衰減比為4:1，記錄此時(shí)的比例系數(shù)Kp為Ks，兩個(gè)峰值之間的時(shí)間周期為周期Ts。
假設(shè)：響應(yīng)曲線(xiàn)峰值衰減比為4:1時(shí)的比例系數(shù) Ks= 4.74;
兩個(gè)峰值之間的時(shí)間周期 Ts = 21.9967；
則，按照上面4:1衰減法控制表中的計(jì)算公式可得：
P控制：比例系數(shù) Kp = Ks = 4.74；
PI控制：比例系數(shù) Kp = Ks/1.2 = 3.95；
積分時(shí)間常數(shù) Ti = 0.5 * Ts = 10.9984;
積分系數(shù) Ki = Kp / Ti = 3.95 / 10.9984 = 0.3591;
P（比例）控制
首先，把控制器設(shè)置成純比例控制，即令積分系數(shù)Ki和微分系數(shù)Kd為零，在北太真元建立模型，形成比例控制系統(tǒng)，結(jié)構(gòu)如下圖所示；
設(shè)置參數(shù)：
仿真時(shí)長(zhǎng)：30s；步長(zhǎng)0.01s；求解器：ode4
得到的仿真結(jié)果，如下圖所示：
PI（比例積分）控制
首先，在純比例控制系統(tǒng)的基礎(chǔ)上增加積分系數(shù)Ki，令微分系數(shù)Kd為零，在北太真元建立模型，形成比例控制系統(tǒng)，結(jié)構(gòu)如下圖所示；
設(shè)置參數(shù)：
仿真時(shí)長(zhǎng)：30s；步長(zhǎng)0.01s；求解器：ode4
得到的仿真結(jié)果，如下圖所示：
公眾號(hào)
版權(quán)所有 ©北太振寰（重慶）科技有限公司 | 渝ICP備2022008375號(hào)-4 法律聲明 隱私政策