彈簧質(zhì)量阻尼實(shí)驗(yàn)指導(dǎo)書(shū)

1、質(zhì)量-彈簧-阻尼系統(tǒng)實(shí)驗(yàn)教學(xué)指導(dǎo)書(shū)Electro mechanic alPlantSyttem InterlaceSoftware("Executive Program*)Input i Output ElectronicsDSP Btittud Controller f Data Acquisition HoardRpsl-tirrw ControJIer A kO北京理工大學(xué)機(jī)械與車(chē)輛學(xué)院2016.3實(shí)驗(yàn)一：?jiǎn)巫杂啥认到y(tǒng)數(shù)學(xué)建模及仿真1實(shí)驗(yàn)?zāi)康?1) 熟悉單自由度質(zhì)量-彈簧-阻尼系統(tǒng)并進(jìn)行數(shù)學(xué)建模;(2) 了解MATLAB軟件編程，學(xué)習(xí)編寫(xiě)系統(tǒng)的仿真代碼;(3) 進(jìn)行單自由度系統(tǒng)

2、的仿真動(dòng)態(tài)響應(yīng)分析。2實(shí)驗(yàn)原理n亠k單自由度質(zhì)量-彈簧-阻尼系統(tǒng)，如上圖所示。由一個(gè)質(zhì)量為 m的滑塊、一個(gè)剛度系數(shù)為k的彈簧和一個(gè)阻尼系數(shù)為c的阻尼器組成。系統(tǒng)輸入：作用在滑塊 上的力f(t)。系統(tǒng)輸出：滑塊的位移x(t)建立力學(xué)平衡方程：? ?mx cx kx f變化為二階系統(tǒng)標(biāo)準(zhǔn)形式：? ?2 Tx 2 x x m其中：3是固有頻率，Z是阻尼比。2m 2ikm2.1欠阻尼(Z1)情況下，輸入f(t)和非零初始狀態(tài)的響應(yīng)：x(t)2 e ")sin( 1 (t )dom Ji2?arctan(j2)2 sin( V1t)2.2欠阻尼(Z1)情況下，輸入f(t)=fo*COS( g*

3、t)和非零初始狀態(tài)的的響應(yīng):f 0m (20Y(223 0f0)2 2 丿0k(2行(2)2訂1x(0)2 0f0(212k J2( 2x(t)x(0)x(0)e t cos( _1輸出振幅和輸入振幅的比值：A =m/2cos( 0t arctan()02t)22 2 2 專(zhuān)亍誌宀(Lt)12 譲(2)2 03動(dòng)力學(xué)仿真根據(jù)數(shù)學(xué)模型，使用龍格庫(kù)塔方法ODE45求解，任意輸入下響應(yīng)結(jié)果。仿真代碼見(jiàn)附件4實(shí)驗(yàn)4.1固有頻率和阻尼實(shí)驗(yàn)(1 )將實(shí)驗(yàn)臺(tái)設(shè)置為單自由度質(zhì)量-彈簧-阻尼系統(tǒng)。(2) 關(guān)閉電控箱開(kāi)關(guān)。點(diǎn)擊setup菜單，選擇Control Algorithm，設(shè)置選擇(3) 點(diǎn)擊Comma

4、nd 菜單，選擇Tectory ,選取step,進(jìn)入set-up,選取OpenLoop Step 設(shè)置(O)counts, dwell time=3000ms,(1)rep, 然后 OK。此步是為了 使控制器得到一段時(shí)間的數(shù)據(jù)，并不會(huì)驅(qū)動(dòng)電機(jī)運(yùn)動(dòng)。(4) 點(diǎn)擊Data菜單，選擇Data Acquisition,設(shè)置選取Encoder#1，然后OK離開(kāi)；從Utility菜單中選擇Zero Position使編碼器歸零。(5 )從Comma nd菜單中選擇Execute，用手將質(zhì)量塊1移動(dòng)到2.5cm左右的 位置(注意不要使質(zhì)量塊碰觸移動(dòng)限位開(kāi)關(guān))，點(diǎn)擊 Run,大約1秒后，放開(kāi)手 使其自由震蕩，

5、在數(shù)據(jù)上傳后點(diǎn)擊 OK。(6) 點(diǎn)擊 Plotting 菜單，選擇 Setup Plot，選取 Encoder #1 Position;然后點(diǎn)擊Plotting菜單，選擇Plot Data，則將顯示質(zhì)量塊1的自由振動(dòng)響應(yīng)曲線(xiàn)。(7) 在得到的自由振動(dòng)響應(yīng)曲線(xiàn)圖上，選擇n個(gè)連續(xù)的振幅明顯的振動(dòng)周期，計(jì)算出這段振動(dòng)的時(shí)間t，由n/t即可得到系統(tǒng)的頻率，將Hz轉(zhuǎn)化為rad/sec即為 系統(tǒng)的振動(dòng)頻率(8) 在自由振動(dòng)響應(yīng)曲線(xiàn)圖上，測(cè)量步驟 7選取時(shí)間段內(nèi)初始振動(dòng)周期的振幅 X0以及末尾振動(dòng)周期的振幅Xn。由對(duì)數(shù)衰減規(guī)律即可求得系統(tǒng)阻尼比。(9) 實(shí)驗(yàn)數(shù)據(jù)記錄序號(hào)第1次實(shí)驗(yàn)第2次實(shí)驗(yàn)第3次實(shí)驗(yàn)實(shí)驗(yàn)測(cè)試

6、頻率實(shí)驗(yàn)測(cè)試阻尼比滑塊質(zhì)量m彈簧剛度k阻尼系數(shù)c頻率理論值阻尼比理論值頻率估計(jì)誤差阻尼比估計(jì)誤差(10)在仿真代碼基礎(chǔ)上，計(jì)算出實(shí)驗(yàn)結(jié)果對(duì)應(yīng)的理論結(jié)果。對(duì)比分析理論和實(shí)驗(yàn)結(jié)果的差異。完成實(shí)驗(yàn)報(bào)告。4.2幅頻特性實(shí)驗(yàn)(1) 點(diǎn)擊 Comma nd 菜單，選擇 Trajectory，選取 Sinu scidal ，進(jìn)入 set-up,選取 Open LoopStep 設(shè)置Amplitude(0.5V), Frequency(2Hz),Repetition(8)，然后OK。(2) 從Utility菜單中選擇Zero Positio n 使編碼器歸零。從Comma nd菜單中選 擇Execute，點(diǎn)擊

7、Run,在數(shù)據(jù)上傳后點(diǎn)擊OK。(3) 然后點(diǎn)擊Plotting菜單，選擇Plot Data，則將顯示滑塊的受迫振動(dòng)響應(yīng)曲 線(xiàn)。在響應(yīng)曲線(xiàn)圖上，測(cè)量出振動(dòng)振幅，計(jì)算出振動(dòng)的頻率并于輸入的正弦曲線(xiàn) 頻率比較。(4) 根據(jù)實(shí)驗(yàn)情況，改變輸入的正弦曲線(xiàn)頻率的大小，重復(fù)上述，紀(jì)錄實(shí)驗(yàn)數(shù) 據(jù)。輸入頻率滑塊實(shí)驗(yàn)幅值滑塊仿真幅值0.1Hz0.2Hz(5) 在仿真代碼基礎(chǔ)上，實(shí)現(xiàn)正弦激勵(lì)代碼，計(jì)算出實(shí)驗(yàn)結(jié)果對(duì)應(yīng)的理論結(jié)果對(duì)比分析理論和實(shí)驗(yàn)結(jié)果的差異。完成實(shí)驗(yàn)報(bào)告實(shí)驗(yàn)二：雙自由度系統(tǒng)數(shù)學(xué)建模及仿真1實(shí)驗(yàn)?zāi)康模?）熟悉雙自由度質(zhì)量-彈簧-阻尼系統(tǒng)并進(jìn)行數(shù)學(xué)建模;（2）了解MATLAB軟件編程，學(xué)習(xí)編寫(xiě)系統(tǒng)的仿真代碼

8、;（3）進(jìn)行雙自由度系統(tǒng)的仿真動(dòng)態(tài)響應(yīng)分析。2實(shí)驗(yàn)原理2.1數(shù)學(xué)建模雙自由度質(zhì)量-彈簧-阻尼系統(tǒng)，如上圖所示。由兩個(gè)質(zhì)量為 口!和m2的滑塊、 兩個(gè)剛度系數(shù)為kl和k2的彈簧和兩個(gè)阻尼系數(shù)為C1和C2的阻尼器組成。系統(tǒng)輸 入：作用在滑塊上的力f（t）。系統(tǒng)輸出：滑塊的位移X1（t）和X2（t）。建立力學(xué)平衡方程：? ?mi X1q X1? ?m2 X2C2 X2?C| X2?G X2k1X1k1X2?C X1 k1 x2f(t)k2x2k1x1 02.2固有頻率將動(dòng)力學(xué)方程寫(xiě)成矩陣形式：m1?0X1ciC1?X1k1k1X110 f(t)x200?m2X2C1 gC2)?X2k1 (k1k2)

9、得到系統(tǒng)的質(zhì)量矩陣M和剛度矩陣K 解行列式可得固有頻率方程：|K 2M | 0可計(jì)算出固有頻率方程：2kikik2/kik2、24kfn1,2 0.5m (一)mim2m1 m2mim2兩個(gè)振動(dòng)模態(tài)，兩個(gè)固有頻率：高模和低模。2.3解耦通過(guò)數(shù)學(xué)變換將微分方程變化為以下形式：? ?' 2yi 2 i i yii yii f(t)? ?2y2 2 2 2 y22 y22 f (t)注意：yi和y2不是滑塊的位移?；瑝K的位移Xi(t)和X2(t)是yi和y2的函數(shù)。3動(dòng)力學(xué)仿真根據(jù)數(shù)學(xué)模型，使用龍格庫(kù)塔方法ODE45求解，任意輸入下響應(yīng)結(jié)果。仿真代碼見(jiàn)附件4實(shí)驗(yàn)4.i固有頻率分析(i )將實(shí)

10、驗(yàn)臺(tái)設(shè)置為雙自由度質(zhì)量-彈簧-阻尼系統(tǒng)，第一個(gè)滑塊沒(méi)有阻尼器可 以不接，認(rèn)為第一個(gè)阻尼為零。(2 )閉合控制器開(kāi)關(guān)，點(diǎn)擊setup菜單，選擇Control Algorithm，設(shè)置選擇Con ti nuous Time Co ntrol，Ts=0.0042,然后 OK。點(diǎn)擊 Comma nd 菜單，選擇Trajectory，選取 step,進(jìn)入 set-up, 選取 Open Loop Step 設(shè)置(0)counts, dwell time=3000ms, (i)rep, 然后OK。此步是為了使控制器得到一段時(shí)間的 數(shù)據(jù)，并不會(huì)驅(qū)動(dòng)電機(jī)運(yùn)動(dòng)。2)點(diǎn)擊 Data 菜單，選擇 Data Acqu

11、isition, 設(shè)置分別選取 Encoder#1Encoder#2 ，然后 OK 離開(kāi)；從 Utility 菜單中選擇 Zero Position 使編碼器歸零。(4 )從Comma nd菜單中選擇Execute，用手將質(zhì)量塊1移動(dòng)到2.5cm左右的 位置(注意不要使質(zhì)量塊碰觸移動(dòng)限位開(kāi)關(guān))，點(diǎn)擊 Run, 大約 1 秒后，放開(kāi)手 使其自由振蕩，在數(shù)據(jù)上傳后點(diǎn)擊 OK 。(5) 點(diǎn)擊 Plotting 菜單，選擇 Setup Plot，分別選取 Encoder #1 Position , Encoder #2Position ；然后點(diǎn)擊 Plotting 菜單，選擇 Plot Data ，則

12、將顯示質(zhì)量塊 1 ， 2的自 由振動(dòng)響應(yīng)曲線(xiàn)。( 6 )實(shí)驗(yàn)數(shù)據(jù)紀(jì)錄：實(shí)驗(yàn)條件：滑塊質(zhì)量m1和m2，彈簧剛度k1和k2，阻尼系數(shù)&和c?。 實(shí)驗(yàn)數(shù)據(jù)：時(shí)間 -滑塊1位移數(shù)據(jù)；時(shí)間 -滑塊2位移數(shù)據(jù)。問(wèn)題1 ：兩個(gè)滑塊位移的頻率測(cè)量值是高模和低模頻率么？ 問(wèn)題2：實(shí)際的機(jī)械系統(tǒng)是多自由度的，如何通過(guò)實(shí)驗(yàn)法測(cè)試系統(tǒng)固有頻率？( 7)實(shí)驗(yàn)報(bào)告。關(guān)鍵點(diǎn)是理論和實(shí)驗(yàn)結(jié)果對(duì)比分析。4.2幅頻特性實(shí)驗(yàn)( 1 )點(diǎn)擊 Command 菜單，選擇 Trajectory ，選取 Sinuscidal, 進(jìn)入 set-up, 選 取 Ope n Loop Step 設(shè)置(200cou nts)Amplitu

13、de, Frequency(2Hz),Repetition ( 8 )，然后 OK 。(2)從Utility菜單中選擇Zero Positio n 使編碼器歸零。從Comma nd菜單中選 擇Execute，點(diǎn)擊Run,在數(shù)據(jù)上傳后點(diǎn)擊OK。3)然后點(diǎn)擊 Plotting 菜單，選擇 Plot Data ，則將顯示質(zhì)量塊 1 ， 2的受迫振 動(dòng)響應(yīng)曲線(xiàn)。在響應(yīng)曲線(xiàn)圖上，即可測(cè)量出振動(dòng)振幅。問(wèn)題1:單自由度和雙自由度系統(tǒng)的幅頻特性有何差異?問(wèn)題2 :高模貢獻(xiàn)分析。實(shí)驗(yàn)三：PID控制1實(shí)驗(yàn)?zāi)康模?）學(xué)習(xí)PID閉環(huán)控制結(jié)構(gòu)和系統(tǒng)閉環(huán)傳遞函數(shù)計(jì)算;（2）PID控制器參數(shù)設(shè)計(jì)；（3 ）控制性能分析。2實(shí)

14、驗(yàn)原理Refarwco Inpm (Eju. Input Trajecloryj上圖給出閉環(huán)控制系統(tǒng)原理框圖。單自由度質(zhì)量 -彈簧-阻尼系統(tǒng)結(jié)構(gòu)下，斷開(kāi)彈簧和阻尼，僅僅保留滑塊質(zhì)量m。?mx f電控箱可以看做比例增益k hw。其中：u是控制器輸出。PID控制：?kpe ki edt kd e u其中：e是比較器輸出，參考輸入與實(shí)際輸出的偏差值。根據(jù)全部上式，可得閉環(huán)結(jié)構(gòu)微分方程:?m x? ?khwkd x khwkp x 氐? ?khwkd r khwkp rkhwkir對(duì)應(yīng)的傳遞函數(shù):x(s) r(s)3 ms(khwkd )s2 (khwkp)s (khwK)(khwkd)s2 (khw

15、kp)s (khwki)3 PID設(shè)計(jì)PID控制器中設(shè)置積分因子ki為零，則為PD控制。傳遞函數(shù)變?yōu)?x(s)khwkdKiwkpsmr(s)m2khwkdkhwkpssm閉環(huán)特征方程是分母:2mkhwkdgj=4 Hzk hw kpkhw k dZ0Z0.2Z1Z2.021” khwkp設(shè)計(jì)頻率 併4Hz，三種阻尼(欠阻尼Z0.2，臨界Z1.0，過(guò)阻尼Z2.0 )的控制器設(shè)計(jì)。4實(shí)驗(yàn) 4.1頻率(1 )在控制器斷開(kāi)的情況下，拆除與質(zhì)量塊1連接的彈簧，使其余元件遠(yuǎn)離質(zhì)量塊1的運(yùn)動(dòng)范圍，為其安裝4個(gè)500g的銅塊，加上小車(chē)本身的質(zhì)量，標(biāo)定總質(zhì)量 m=2.6kg(2)實(shí)驗(yàn)標(biāo)定khw值：根據(jù)估計(jì)出的

16、khw值，設(shè)置控制器ki=0和kd=0，調(diào)整kp來(lái) 估計(jì)系統(tǒng)系統(tǒng)頻率w=4Hz。注意：kp不能大于0.08。( 3)閉合控制器開(kāi)關(guān)；點(diǎn)擊 Data 菜單，選擇 Data Acquisition, 設(shè)置選取En coder#1 和 Comma nded Positi on in formatio n; 點(diǎn)擊 Comma nd 菜單，選擇Tectory，選取 step，設(shè)置(0)counts,dwell time=3000ms,rep 。( 4)點(diǎn)擊 setup 菜單，選擇 Control Algorithm ，設(shè)置選擇 Continuous TimeControl ，Ts=0.0042,選取PI

17、D，進(jìn)入Setup Algorithm ，輸入 kp的值(ki=0和 kd=0)(輸入的值不能大于kp=0.08 )，然后OK。移動(dòng)質(zhì)量塊1到-0.5cm的位 置(規(guī)定，朝電機(jī)方向?yàn)樨?fù))選擇Implement Algorithm ，然后Ok。注意：從此步開(kāi)始的每一步， 要進(jìn)行下一步之前， 都要與運(yùn)動(dòng)裝置保持一定的安全距離;選擇 Implement Algorithm 后，控制器將會(huì)立即加載，若出現(xiàn)不穩(wěn)定 的或者很大的控制信號(hào)時(shí)， 運(yùn)動(dòng)裝置可能反應(yīng)很劇烈; 若加載后， 系統(tǒng)看上去穩(wěn) 定，要先用一輕質(zhì)不尖銳的物體輕輕碰觸質(zhì)量塊以驗(yàn)證其穩(wěn)定性。( 5)點(diǎn)擊 Command 菜單，進(jìn)入 Execute

18、 ，用手將質(zhì)量塊移動(dòng)到 2cm 左右的位 置，點(diǎn)擊 Run, 移動(dòng)質(zhì)量塊大約到 3cm 位置，然后釋放(不要拿著質(zhì)量塊多于一 秒，以免電機(jī)過(guò)熱而斷開(kāi)控制)。(6) 點(diǎn)擊 Plotting 菜單，選擇 Setup Plot，選取Encoder #1 ;然后點(diǎn)擊 Plotting菜單，選擇 Plot Data ，則將顯示質(zhì)量塊 1的時(shí)間響應(yīng)曲線(xiàn)。試想一下，若將比例增益系數(shù)kp增加一倍，則系統(tǒng)的響應(yīng)頻率將有什么變化？4.2 阻尼(1)確定kd的值(不能大于0.04 )，使得khwkd=50N/(m/s)，重復(fù)步驟4，除 了輸入kd的值和ki=0和kp=O(2) 先用尺子檢查系統(tǒng)的穩(wěn)定性，然后用手來(lái)回

19、的移動(dòng)質(zhì)量塊來(lái)感受系數(shù)kd帶 來(lái)的粘性阻尼的影響 (注意不要極度的迫使質(zhì)量塊運(yùn)動(dòng)， 以免電機(jī)過(guò)熱而斷開(kāi)控 制)。(3) 增大的kd值(kd<0.04 )，重復(fù)以上步驟，看能否感受到阻尼的增大？4.3 位置控制( 1)控制電機(jī)驅(qū)動(dòng)滑塊 1 移動(dòng)到某一個(gè)特定的位置，獲取時(shí)間 -位移數(shù)據(jù)后， 計(jì)算閉環(huán)控制特性，包括：超調(diào)量、上升時(shí)間、調(diào)節(jié)時(shí)間、穩(wěn)態(tài)誤差。( 2)使用仿真方法計(jì)算得到同樣條件下的動(dòng)態(tài)響應(yīng)結(jié)果，對(duì)比分析理論與實(shí)驗(yàn) 結(jié)果的差異。4.4 正弦激勵(lì)(1)使用正弦輸入驅(qū)動(dòng)滑塊作正弦運(yùn)動(dòng)， 獲取時(shí)間 -位移數(shù)據(jù)后， 計(jì)算閉環(huán)控制 系統(tǒng)幅頻特性。( 2)使用仿真方法計(jì)算得到同樣條件下的幅頻響應(yīng)

20、結(jié)果，對(duì)比分析理論與實(shí)驗(yàn) 結(jié)果的差異。4.5 柔性結(jié)構(gòu)控制( 1)實(shí)驗(yàn)臺(tái)只有質(zhì)量塊，為剛性結(jié)構(gòu)機(jī)械系統(tǒng)。如果安裝上一根彈簧，則組成 簡(jiǎn)單的柔性結(jié)構(gòu)系統(tǒng) (具有柔性環(huán)節(jié)) 。此時(shí)機(jī)械系統(tǒng)由于能量交互可能產(chǎn)生振 動(dòng)。如何設(shè)計(jì)控制器來(lái)動(dòng)態(tài)控制機(jī)械振動(dòng)。實(shí)驗(yàn)四：擾動(dòng)抑制1實(shí)驗(yàn)?zāi)康?1) 外擾動(dòng)設(shè)置；(2) PID控制器設(shè)計(jì);(3 )控制性能分析。2實(shí)驗(yàn)原理r上圖給出閉環(huán)控制系統(tǒng)原理框圖。單自由度質(zhì)量 -彈簧-阻尼系統(tǒng)結(jié)構(gòu)下，斷開(kāi)彈簧和阻尼，僅僅保留滑塊質(zhì)量 m。位移輸出受到外擾動(dòng)d干擾?m x電控箱可以看做比例增益k hw其中：u是控制器輸出。PID控制:kpe ki edt kd e u其中：e是比

21、較器輸出，參考輸入與實(shí)際輸出的偏差值。根據(jù)全部上式，可得閉環(huán)結(jié)構(gòu)微分方程:? ? ?m xkhwkd x khwkpx khwx? ? ?khwkd rkhwkp r khwkird對(duì)應(yīng)的傳遞函數(shù):x(s)(khwkd )s2 (khwkp)s (khwki)ms3 (khwkd)s2 (khwkp)s(khwk)(S)ms3 (khwkd )s2(khwkp)s (khwki)d(s)3 PID設(shè)計(jì)PID控制器中設(shè)置積分因子ki為零，則為PD控制。傳遞函數(shù)變?yōu)?x(s) r(s)khwkd ckhwkps -m m閉環(huán)特征方程是分母:khw kp mkhwKjkhw kd2m 2'k

22、hwkpm附件1：實(shí)驗(yàn)1仿真代碼% exp 1clcclearglobal SK SC SM SF;sampli ng=1/1000;SK=200;SC=0;SM=2.6; State=zeros(1,2);State(1)=1; State(2)=0.0; sss=zeros(1,1);for k=1:fix(5/sampling) t=k*sampling;sss(k,1)=t; %timeSF=0; sss(k,2)=SF;TimeOdeArray,VarOdeArray = ode45( mdlDerivatives1, t t+sampling,State); m,n=size(Ti

23、meOdeArray);TimeAtEndOfArray = TimeOdeArray(m,1); if( abs(TimeAtEndOfArray - (t+sampling) ) >= abs(0.001*sampling) ) warning( 'numerical integration failed' ); break ; endVAR = VarOdeArray(m, : );%o ?2?3 O ? de ?卩? x'i?State=VAR; %3? e ? x'i ? u D? £ ?DD? o ? ' ? yu'd

24、 sss(k,3)=VAR(1);% responsesss(k,4)=VAR(2);% velocityendplot(sss(:,1),sss(:,2),'r',sss(:,1),sss(:,3),'k' ,sss(:,1),sss(:,4),'b');% mdlDerivatives1function dx=mdlDerivatives1(T,x) global SK SC SM SF;% y(1)=x;% y(2)=xd; dx=zeros(2,1);%sloshing dynamicsdx(1)=x(2); dx(2)=SF/SM-SK

25、/SM*x(1)-SC/SM*x(2);附件 2: 實(shí)驗(yàn) 2 仿真代碼% exp 2clcclearglobal SK1 SK2 SC1 SC2 SM1 SM2 SF;sampling=1/1000;SK1=200;SK2=200;SC1=0;SC2=0;SM1=2.6;SM2=2.6;State=zeros(1,4);State(1)=1.0;State(3)=0.0;sss=zeros(1,1);for k=1:fix(5/sampling) t=k*sampling; sss(k,1)=t; %timeSF=0; sss(k,2)=SF;TimeOdeArray,VarOdeArray

26、= ode45( mdlDerivatives2, t t+sampling,State); m,n=size(TimeOdeArray);TimeAtEndOfArray = TimeOdeArray(m,1);if( abs(TimeAtEndOfArray - (t+sampling) ) >= abs(0.001*sampling) ) warning( 'numerical integration failed' ); break ; endVAR = VarOdeArray(m, : );%o ?2?3 O ? de ?卩? x'i?State=VAR; %3?e? x'i ? UD?£ ?DD? o ? ' ?小'usss(k,3)=VAR(1);% response1sss(k,4)=VAR(2);% velocity1sss(k,5)=VAR(3);% response2sss(k,6)=VAR(4);% velocity2endplot(sss(:,1),sss(:,2),'r' ,sss