攔截飛彈的發射時間窗口決定

國立交通大學

電控工程研究所

碩 士 論 文

攔截飛彈的發射時間窗口決定

Launch Time Window Determination for Missile Interceptor

研 究 生 ： 林彥伯

指 導 教 授 ： 林清安 教授

攔截飛彈的發射時間窗口決定

Launch Time Window Determination for Missile Interceptor

研 究 生：林彥伯 Student: Yen-Po Lin

指導教授：林清安 Advisor: Prof. Ching-An Lin

國 立 交 通 大 學

電控工程研究所

碩 士 論 文

A Thesis

Submitted to Institute of Electrical and Control Engineering College of Electrical and Computer Engineering

National Chiao Tung University In Partial Fulfillment of the Requirements

For the Degree of Master In

Electrical and Control Engineering October 2010

Hsinchu, Taiwan, Republic of China

攔截飛彈的發射時間窗口決定

學生：林彥伯 指導教授：林清安 教授

國立交通大學電控工程研究所

中文摘要

本論文探討攔截飛彈的發射時間窗口(Launch Time Window)該如

何適當選擇，如何描述飛彈及目標的位置與速度函數，以及分析攔截

成功的必要條件並建立最佳化問題，透過發展處罰函數(PF)演算法及

連續二次規劃(SQP)演算法，來求出發射前等待的延遲時間最大值及

最小值，即發射時間窗口。最後提供三種不同來襲目標來驗證所求出

的發射時間窗口是否符合實際攔截情形。透過模擬的數值結果顯示針

對不同目標所求出的發射時間窗口均有不錯的準確度且亦相當符合

實際攔截情形。

Launch Time Window Determination for Missile Interceptor

Student: Yen-Po Lin Advisor: Prof. Ching-An Lin

Institute of Electrical and Control Engineering

National Chiao Tung University

Abstract

When an attacking missile is detected, the missile defense system

must decide when to launch the missile interceptor. The choice of a

proper launch time is important because given the physical constraints on

the interceptor, a poor choice of launch time alone may cause the

interception to fail. This thesis discusses the determination of launch time

window for a missile interceptor, which is formulated as constrained

optimization problems.

The earliest and the latest allowable launch times are determined by

numerically solve the corresponding minimization and maximization

problems. Three incoming threats are considered. Numerical results show

good agreement with results obtained via six degree of freedom

誌 謝

本論文可以順利完成，首先要感謝我的指導教授林清安博士在研

究上提供我許多寶貴的經驗及建議，使我得以能在研究上順利進行。

不僅僅是在課業上的幫助，在待人處世及面對問題該有正確、積極的

態度也是我值得學習的榜樣。接著感謝口詴委員鄧清政教授與陳科祥

博士仔細地閱讀我的論文並給予我適當的指教。

此外也感謝所有在求學過程中幫助我的朋友，特別是606實驗室

的吳建賢與林世修學長，同窗邱泰偉、蕭惟庭及學弟卓國展在日常生

活上的照顧及學習上的討論，陪伴我渡過這兩年充實且美好的歲月。

最後也要特別感謝我的家人無怨無悔的為我付出，使我衣食無

缺，讓我更有學習的動力。謝謝我的女友佩蓉默默的支持我，鼓勵我。

謝謝交通大學提供我如此好的學習環境，讓我能在專業領域更上層

樓。

目錄

中文摘要... i 英文摘要... ii 誌 謝... iii 目錄... iv 表目錄... vii 圖目錄... ix 第一章、緒論... 1 第二章、運動方程式... 3 2.1 座標系統... 3 2.2 飛彈六自由度運動方程式... 4 2.2.1 平移運動方程式 ... 4 2.2.2 轉動方程式 ... 5 2.2.3 尤拉角(Euler angle) ... 5 2.2.4 尤拉方程式(Euler equation) ... 7 2.3 目標運動方程式... 7 第三章、飛彈及目標的速度及位置函數... 9 3.1 各階段時間的定義... 9 3.2 變數的定義... 10 3.3 目標的速度及位置函數... 11 3.4 飛彈的速度及位置函數... 11 3.5 第二節火箭加速度分析... 13 第四章、最佳化問題的建構... 15 4.1 攔截的必要條件... 15 4.2 最佳化問題的建立... 18

4.3 延遲時間最大值與最小值... 19 第五章、最佳化演算法的發展... 21 5.1 處罰函數(PF)演算法 ... 21 5.2 PF 範例 ... 23 5.2.1 PF 範例最小值問題 ... 23 5.2.2 PF 範例最大值問題 ... 25 5.3 連續二次規劃(SQP)演算法... 27 5.4 SQP 範例 ... 29 5.4.1 SQP 範例最小值問題 ... 29 5.4.2 SQP 範例最大值問題 ... 30 第六章、數值計算及模擬... 33 6.1 飛彈的推力及物性... 33 6.2 模擬初始條件... 34 6.3 模擬結果... 39 6.3.1 Simulink 模擬目標一來襲 ... 40 6.3.1.1 執行 PF 演算法求解目標一 ... 44 6.3.1.2 執行 SQP 演算法求解目標一 ... 46 6.3.2 Simulink 模擬目標二來襲 ... 47 6.3.2.1 執行 PF 演算法求解目標二 ... 50 6.3.2.2 執行 SQP 演算法求解目標二 ... 52 6.3.3 Simulink 模擬目標三來襲 ... 53 6.3.3.1 執行 PF 演算法求解目標三 ... 55 6.3.3.2 執行 SQP 演算法求解目標三 ... 57 6.4 延遲時間評估與驗證... 57 6.4.1 當目標一來襲時 ... 59

6.4.3 當目標三來襲時 ... 64 6.5 不同h 參數的數值結果 ... 67_r 第七章、結論... 71 參考資料... 72 附錄 1 主程式碼... 73 附錄 2 限制式程式碼... 74 附錄 3 PF 程式碼 ... 76 附錄 4 SQP 程式碼 ... 86 附錄 5 SQP 求解不同目標的 LNS ... 91

表目錄

表 2.1：慣性座標轉換至體座標的旋轉順序... 6 表 4.1：選擇不同延遲時間對攔截的影響... 20 表 5.1：PF 求解範例最小值相關數據 ... 24 表 5.2：PF 求解範例最大值相關數據 ... 26 表 5.3：SQP 求解範例最小值相關數據 ... 30 表 5.4：SQP 求解範例最大值相關數據 ... 31 表 5.5：最小值問題方法比較... 31 表 5.6：最大值問題方法比較... 31 表 6.1：比較探空五號與本論文火箭的推力與作用時間... 33 表 6.2：比較探空五號與本論文火箭的燃料質量... 34 表 6.3：飛彈燃料質量分配狀況... 34 表 6.4：目標資料庫... 34 表 6.5：模擬目標一在不同延遲時間的攔截情形... 41 表 6.6：PF 求解目標一延遲時間最大值 ... 45 表 6.7：PF 求解目標一延遲時間最小值 ... 46 表 6.8：SQP 求解目標一延遲時間最大值 ... 46 表 6.9：SQP 求解目標一延遲時間最小值 ... 47 表 6.10：模擬目標二在不同延遲時間的攔截情形... 48 表 6.11：PF 求解目標二延遲時間最大值 ... 51 表 6.12：PF 求解目標二延遲時間最小值 ... 52 表 6.13：SQP 求解目標二延遲時間最大值 ... 52 表 6.14：SQP 求解目標二延遲時間最小值 ... 53 表 6.15：模擬目標三在不同延遲時間的攔截情形... 54 表 6.16：PF 求解目標三延遲時間最大值 ... 56

表 6.17：SQP 求解目標三延遲時間最大值 ... 57 表 6.18：三種不同來襲目標的延遲時間窗口與可靠攔截時間... 58 表 6.19：目標一求解可靠攔截時間代入模擬結果... 59 表 6.20：目標二求解可靠攔截時間代入模擬結果... 62 表 6.21：目標三求解可靠攔截時間代入模擬結果... 64 表 6.22：當HR=5 求解延遲時間最小值 ... 67 表 6.23：當HR=10 求解延遲時間最小值 ... 68 表 6.24：不同高度差產生的延遲時間誤差及可靠攔截時間... 69

圖目錄

圖 2.1：NED 座標示意圖 ... 3 圖 2.2：尤拉角示意圖... 6 圖 3.1：各時間參數的相對關係... 10 圖 4.1：目標與飛彈對座標原點的距離示意圖... 16 圖 5.1：PF 範例最小值驗証 ... 24 圖 5.2：PF 範例最大值驗証 ... 26 圖 6.1：飛行過程中的質量與推力隨時間變化示意圖... 33 圖 6.2：不同目標的飛行軌跡圖(NED) ... 35 圖 6.3：不同目標的飛行高度對時間圖... 36 圖 6.4：飛彈與目標一的飛行軌跡圖(NED) ... 42 圖 6.5：飛彈與目標一的飛行高度對時間圖... 42 圖 6.6：目標一不同延遲時間對應的零力誤失距離圖... 43 圖 6.7：目標一不同延遲時間對應的高度關係圖... 44 圖 6.8：目標一延遲時間示意圖... 45 圖 6.9：飛彈與目標二的飛行軌跡圖(NED) ... 49 圖 6.10：目標二不同延遲時間對應的零力誤失距離圖... 50 圖 6.11：目標二延遲時間示意圖 ... 51 圖 6.12：飛彈與目標三的飛行軌跡圖(NED) ... 55 圖 6.13：飛彈與目標三的飛行高度對時間圖... 55 圖 6.14：目標三延遲時間示意圖... 56 圖 6.15：驗證飛彈與目標一的飛行軌跡圖(NED) ... 60 圖 6.16：驗證飛彈與目標一的飛行高度對時間圖... 61 圖 6.17：驗證飛彈與目標一對座標原點的距離圖... 61 圖 6.18：驗證飛彈與目標二的飛行軌跡圖(NED) ... 63

圖 6.19：驗證飛彈與目標二的飛行高度對時間圖... 63 圖 6.20：驗證飛彈與目標二對座標原點的距離圖... 64 圖 6.21：驗證飛彈與目標三的飛行軌跡圖(NED) ... 65 圖 6.22：驗證飛彈與目標三的飛行高度對時間圖... 66 圖 6.23：驗證飛彈與目標三對座標原點的距離圖... 66 圖 6.24：不同高度差與實際延遲時間所對應的零力誤失距離圖... 69

第一章

緒論

早在 1995 及 1996 年台海危機時，中國假藉軍事演習之名發射飛彈至台灣 附近海域並完成五階段「攻台模式」演練，包括發射彈道飛彈摧毀機場、制空演 習、制海演習、登陸演習及城鎮山地戰，如此大膽的挑釁行為，使得國人日漸感 受到中國飛彈的真實威脅。特別是近年來，中國積極部署對台飛彈，更多次阻止 美國對台軍售，不僅在沿海地區部署超過一千多枚飛彈對準台灣，且隨著軍事預 算年年增加，更在內陸的雲南及廣西柳州部署好幾百顆高精準度且射程不等的彈 道飛彈，可以從不同的角度與距離瞄準台灣，其射程最遠超過 1500 公里。面對 軍事武力如此強大的中國，我國更應積極發展攔截飛彈防衛系統，特別是面對射 程不同的目標，該有不同的發射時間點，因此本論文主要探討當敵方朝我方發射 飛彈來攻擊某一目標時，我方應如何選擇最佳發射時間來攔截目標飛彈。 當地面雷達站偵測到有一目標從遠方發射至我方附近領土時，在典型的飛 彈防禦系統中，雷達偵測到來襲目標的距離通常遠大於防禦飛彈的射程，例如飛 彈的有效射程在 200 至 300 公里而雷達偵測的能力超過 1000 公里。在這種狀況 下，偵測到目標如果立即發射防禦飛彈，攔截必定失敗。所以我們需要有一個延 遲時間t 可以發射攔截飛彈使得在第二節中途導引[7]結束後產生足夠小的零力_d

誤失距離(Zero Effort Miss，ZEM)，讓第三節截殺彈頭在與目標相對距離夠近時 鎖定並開始追蹤目標，以致命中目標；而且攔截的高度通常也有一定的下限。因 此飛彈發射的時間必頇適當地選擇才能確保攔截成功。

如何在有限的時間範圍內進行攔截即飛彈可以發射的時間窗口(Launch Time Window)，將是本論文需要深入探討的問題。在探討過程中，為了計算上的 方便與一致性，位置的單位皆使用公里(km)，速度的單位皆使用公里/每秒 (km/s)，加速度的單位皆使用公里/每秒平方(km/s2)。 本論文討論的重點在於如何描述飛彈及目標的速度及位置函數，以及在滑 行時間固定的狀況下如何建立最佳化問題來求出發射前等待的延遲時間最大值 及最小值，此上下限定義了飛彈可以發射的時間窗口。之後介紹兩種演算法來針 對建構出的最佳化問題進行求解，分別是處罰函數方法(Penalty Function,PF method)[3][5] 及 連 續 二 次 規 劃 的 方 法 (Sequential Quadratic Programming,SQP method)[1][2][4]，並提供一個簡單的非線性最佳化問題來測詴所提出的 PF 及 SQP 演算法是否具有可行性與比較兩種方法的差異性。最後提供三種不同來襲目 標來說明延遲時間的重要性以及利用 PF 與 SQP 演算法求出的延遲時間窗口及飛 彈飛行時間來驗證是否接近實際攔截的情況。 本論文第二章介紹運動方程式，第三章推導飛彈及目標的速度及位置函 數，第四章建構最佳化問題，第五章介紹最佳化演算法，第六章為數值計算及模 擬，第七章為結論。

第二章

運動方程式

飛彈在空間中的運動可分為平移運動與轉動，各為三個自由度的運動，因此 在整個三維空間中視飛彈為六個自由度的運動。而將目標假設為點質量，不考慮 其姿態變化，故只考慮目標在三維空間中的平移運動。本章先介紹所需要的座標 系統，再接著介紹飛彈六自由度運動方程式及目標運動方程式。

2.1 座標系統

為了方便描述攔截飛彈與目標在三維空間飛行過程中的動態與姿態，我們需 要定義以下座標系統。

地面座標系( North-East-Down coordinate, NED )

使用地面座標是為了方便描述攔截飛彈與目標相對於地表的軌跡。地面座標

是以地表的雷達站當作座標原點，其X 軸指向地表平面雷達站的北方，其_N Y 軸_N

指向地表平面雷達站的東方，其Z 軸滿足右手定則，指向地下，如圖 2.1 所示。 _N

本論文以地面座標系當作慣性座標。

2.2 飛彈六自由度運動方程式

2.2.1 平移運動方程式 平移運動用來描述飛彈質心在慣性座標中的位置變化，平移的運動方程式為 mx mx my my mz mz v a v a v a    (2.1) mx mx my my mz mz r v r v r v    (2.2) 其中



a_mx,a_my,a_mz



、



v_mx,v_my,v_mz



、



r_mx,r_my,r_mz



分別為加速度、速度和位置在慣 性座標中的三個分量。簡單來說，位置的微分為速度，速度的微分為加速度。 飛彈加速度包含飛彈的引擎推力造成的加速度、氣動力造成的加速度以及重 力加速度 1 mx Tx Ax x my NB Ty Ay y mz Tz Az z a F F g a T F F g mass a F F g         _  _   _                    (2.3) 其中 mass 為飛彈質量，TNB為體座標轉換到慣性座標的轉換矩陣(如 2.2.3 節所討 論)，(F_Tx,F_Ty,F_Tz)、(F_Ax,F_Ay,F_Az)分別為引擎推力和氣動力在體座標中的三個分 量，(g g g_x, _y, _z)為重力加速度在慣性座標中的三個分量。

2.2.2 轉動方程式 旋轉運動在體座標中表示，轉動運動方程式可表示如下













x z y x x y x z y y z y x z z M qr p I I I I M rp q I I I I M pq r I I I I          (2.4) 其中(M M_x, _y,M_z)、( ,I I I_{x y}, _z)、( , , )p q r 分別為力矩、轉動慣量、角速度在體座 標中的三個分量。 力矩包含了引擎推力造成的力矩和氣動力造成的力矩 x Tx Ax y Ty Ay z Tz Az M M M M M M M M M        _{ }               (2.5) 其中(M_Tx,M_Ty,M_Tz)為引擎力矩在體座標中的三個分量，(M_Ax,M_Ay,M_Az)為氣動 力矩在體座標中的三個分量。 2.2.3 尤拉角(Euler angle) 慣性座標與體座標之間的關係可以透過一組尤拉角



  , ,



來表示，通常用來描 述一剛體的方位。如圖 2.2 中，首先將慣性座標與體座標的原點重合，X 為體_B 座標的 X 軸，zB為體座標Z 軸在B Y Z 平面上的投影，N- N  、、分別為三個有 順序的旋轉角度。

  XN N Y N Z B x  B z   XN N Y N Z B x  B z 圖 2.2：尤拉角示意圖 在此將定義由慣性座標轉換至體座標的旋轉順序整理成如表 2.1 所示： 順序 1 以Z 軸為轉軸並旋轉_N  角 順序 2 以此時的Y 為轉軸並旋轉_N 角 順序 3 以此時的X 為轉軸並旋轉_B 角 表 2.1：慣性座標轉換至體座標的旋轉順序 經過表 2.1 有順序的三個旋轉後，即可將慣性座標旋轉到體座標的方向，因此兩 者之間的轉換矩陣即為三個旋轉矩陣有順序的相乘而得到。從慣性座標轉換到體 座標的轉換矩陣為

1 0 0 cos 0 sin cos sin 0

0 cos sin 0 1 0 sin cos 0

0 sin cos sin 0 cos 0 0 1

cos cos cos sin sin

cos sin sin cos sin sin sin sin cos cos sin cos

cos cos sin sin sin sin cos

BN T                                                 _{    } _                 

 sin sin cos  cos cos 

           (2.6)

若是要計算從體座標轉換到慣性座標的轉換矩陣則為 ( )T NB BN T  T (2.7) 其中 (TBN)T為TBN的轉置矩陣。 2.2.4 尤拉方程式(Euler equation) 我們假設( , , )p q r 分別為飛彈在滾轉、俯仰及偏航方向的角速度，再搭配前 ㄧ節所定義的尤拉角( , , )   ，透過簡單的推導可以得到尤拉角和角速度之間的 一階微分方程式，稱為尤拉方程式[10]：









tan sin cos

cos sin

sec sin cos

p q r q r q r                   (2.8)

2.3 目標運動方程式

假設目標為點質量，不考慮其姿態變化，故只考慮目標在三維空間中的平移 運動，其運動方程式在慣性座標下可表示為 tx tx ty ty tz tz v a v a v a    (2.9) tx tx ty ty tz tz r v r v r v    (2.10) 其中



a a a_tx, _ty, _tz



、



v v v_tx, _ty, _tz



、



r r r_tx, _ty, _tz



分別為加速度、速度和位置在慣性座標

考慮目標加速度只受到重力和空氣阻力的影響，在慣性座標下可表示為 tz 0 0 tx Dx ty Dy tz Dz a A a A a A g         _   _                   (2.11) 其中(A_Dx,A_Dy,A_Dz)為空氣阻力所產生的加速度在慣性座標下的三個分量，g 為_tz 目標重力加速度在慣性座標下的 Z 分量。因空氣阻力所產生之加速度與目標的 高度，即空氣密度，及馬赫數，即速度，有關。在數值模擬的例子中，我們採用 參考資料[8]有建立的模式。

第三章

飛彈及目標的速度及位置函數

在本章先說明飛彈發射前各階段時間的定義，接著介紹推導飛彈及目標的速 度及位置函數所需要的參數及變數，之後再推導出飛彈及目標的速度及位置函 數。飛彈以固定的仰角及方向角發射，第一節燃畢前的位置及速度可由六自由度 模擬準確地得到，第二節燃燒期間的加速度以簡單的時間一階函數近似，假設飛 彈在其餘的時間的動態僅受重力的影響，則飛彈在終端航程的位置可方便地表示 成時間的二次函數，方便最佳化計算。模擬的結果顯示，簡化的(無導引)飛彈軌 跡與六自由度模擬中途導引的軌跡差異不是很大。

3.1 各階段時間的定義

我們假設飛彈發射的時間為時間 0 秒，即時間座標原點，而雷達偵測到目標 的時間為t_d秒。定義以下六個相關的時間參數，各時間參數的相對關係如圖 3.1 所示。時間的單位是秒(s)。 d t ：雷達偵測到目標至飛彈發射的延遲時間。 1 t ：第一節火箭的燃燒時間。 2 t ：第二節火箭的燃燒時間。 c t ：滑行時間。 h t ：第二節火箭脫離到第三節姿態控制開始的時間。 f t ：預估攔截飛彈擊中目標的時間。

圖 3.1：各時間參數的相對關係

3.2 變數的定義

以下定義描述目標與飛彈運動的變數，距離的單位是公里(km)，時間的單位 是秒(s)，質量的單位是公斤(kg)，力量的單位是牛頓(N)。 0 mx r 、r_my0、rmz0：攔截飛彈初始位置，分別為慣性座標中的三個分量。 mx r 、r_my、r ：攔截飛彈所在位置，分別為慣性座標中的三個分量。 mz tx r 、r_ty、r ：目標所在位置，分別為慣性座標中的三個分量。 _tz 1 V  ：第一節火箭燃燒完畢後的速度，為慣性座標中的三維向量。 1 R  ：第一節火箭燃燒完畢後的位置，為慣性座標中的三維向量。 0 t v ：目標初始速度，為慣性座標中的三維向量。 0 t r ：目標初始位置，為慣性座標中的三維向量。 g：重力加速度，為慣性座標中的三維向量。 2 T ：第二節火箭燃燒時的最大推力。 20 m ：第一節火箭脫離後的質量。 2 m ：第二節火箭的質量燃燒速率。

3.3 目標的速度及位置函數

假設目標的初始位置為r ，初始速度為_t0 v ，只考慮重力的影響，則在任意t0 時間t t_d，目標速度為 0 ( ) ( ) t t d v t v   g t t (3.1) 目標位置為 2 0 0 1 ( ) ( ) ( ) 2 t t t d d r t r v  t t  g t t (3.2) 其中g為重力加速度向量，g(0, 0, 0.0098)km/s2。由(3.2)式可以看出目標的位置 為 t 及t 的二次函數。由於目標的飛行在攔截的過程中均在 30 公里以上的高空，_d 氣動力的影響不大，(3.2)式的模式與實際三自由度模擬的結果差異不大。

3.4 飛彈的速度及位置函數

地面雷達站位於慣性座標原點(0, 0, 0)。假設飛彈發射點的位置為 0 0 0 (r_mx ,r_my ,r_mz )，則高度為rmz0公里。以下描述飛彈在每一個飛行階段的位置及 速度函數。 令第一節火箭燃畢，即時間tt1，飛彈的位置與速度分別為R1與V1。由 於飛彈的推力物性及氣動模式固定，R₁及V₁只和飛彈發射的角度有關。在不 同的發射角度下，這兩個向量可利用六自由度模擬很準確且很有效率的取得。 當時間t t1 tc，即當滑行結束第二節火箭點火前，飛彈的位置與速度分別 為 2 1 1 1 1 ( ) 2 m c c c r t t    R V t  gt (3.3)

1 1 ( ) m c c v t t   V gt (3.4) (3.3)式及(3.4)式的計算只考慮重力忽略了氣動力的影響因此有些誤差。如果t 為c 固定值，則r t_m(₁t_c)及v t_m(₁t_c)可以利用模擬比較準確的取得。大致來說，(3.3) 式及(3.4)式所計算的數值會略高於模擬所得到的數值。 第二節火箭點火之後至燃畢前，即t₁    t_c t t₁ t_c t₂，推力所產生的加速 度因質量遞減而遞增。當中途導引律啟動時，速度的方向隨攔截狀況而改變，因 此加速度的方向並不固定。為了分析方便，也因為實際的加速度方向無從得知， 假設加速度的方向固定並且可近似為線性函數 1 ( ) ( ) m c a t    a t t t b (3.5) 其中 a 與 b 為常數向量，其決定在 3.5 節介紹。因此在第二節點火時間內飛彈的 位置與速度分別為 2 2 3 1 1 1 1 1 1 1 1 1 ( ) ( ) ( ) ( ) ( ) 2 2 6 m c c r t    R V tt  g tt  t t t  a t t t b (3.6) 2 1 1 1 1 1 ( ) ( ) ( ) ( ) 2 m c c v t   V g tt   t t t  a t t t b (3.7) 將t  t1 tc t2代入(3.6)式及(3.7)式，可得時間t  t1 tc t2時飛彈的位置與 速度分別為 2 2 3 2 1 1 2 2 2 2 1 1 1 ( ) ( ) 2 2 6 m f c c r     R V t t  g t t  t  a t b (3.8) 2 2 1 2 2 2 1 ( ) 2 m f c v   V g t t   t a t b (3.9) (3.8)式與(3.9)式分別為飛彈在第二節火箭燃畢時的位置與速度。

當t  t₁ t_c t₂，即截殺彈頭終端飛行階段，其位置與速度分別為 2 2 2 1 2 1 2 1 ( ) ( ) ( ) 2 m m f m f c c r t r v    t t t t  g t  t t t (3.10) 2 1 2 ( ) ( ) m m f c v t v g t  t t t (3.11) 其中r_{m f2} 及vm f2 分別定義於(3.8)式及(3.9)式。

3.5 第二節火箭加速度分析

在真正的攔截飛行中，飛彈中途階段的加速度方向由中途導引律提供，但是 在射控階段，則必頇對加速度做一些簡單的假設以利攔截時間的分析。我們假設 飛彈在滑行完畢中途飛行時，飛行路徑角維持不變且攻角為零。因此在這個階 段，即t₁    t_c t t₁ t_c t₂，飛彈的加速度為 2 1 20 1 2

( )

(

)

m V c

T

a t

L

m

t

t

t

m





  



(3.12) 1 1 1 c V c

V

gt

L

V

gt

 



 

(3.13) 其中T 為第二節燃燒時的最大推力，2 m 為第一節脫離後的質量，常數20 m 為第二2 節火箭的燃燒速率，L 為第一節滑行時間結束後速度的單位向量。 _V1 在滑行時間t 固定的狀況下，在(3.13)式中，_c L 為一常數向量，但是_V1 a t 並_m( ) 非時間 t 的線性函數；為了計算方便，可將(3.12)式右邊作線性化得 2 1 20 ( 1 ) 2 V c T L m   t t t m 

2 1 2 20 1 20 1 1 ( ) V c T L m m t t t m      2 1 2 1 20 20 1 ( ) V c T L m t t t m m    _    _   2 1 2 2 1 1 2 20 20 ( ) V V c T L T m L t t t m m      (3.14) 因此(3.5)式可表示成時間的線性函數為 1 ( ) ( ) m c a t    a t t t b 其中 2 1 20 V T a L m  ， 2 2 1 2 20 V T m b L m 

第四章

最佳化問題的建構

在本章首先說明飛彈能夠攔截成功的必要條件，隨後介紹如何建立一個有效 的非線性最佳化問題來求出發射前的延遲時間及飛彈飛行時間，並透過 Matlab 程式撰寫求出的延遲時間的最大值與最小值，根據此上下限將可以建立發射時間 窗口。最後透過表 4.1 可以初步了解延遲時間的選擇對攔截是否成功的影響。

4.1 攔截的必要條件

一個成功的攔截必頇經過準確的中途導引及終端導引階段，而在發射之初無 法預知中途及終端階段的飛彈加速度。因此如果要探討攔截是否可能，數值上必 頇對飛彈的加速度做一些簡化的假設以利軌跡及速度函數的建立。此即第三章所 建立之結果。 一個成功的攔截的充分條件是在某一個時間tf 0，飛彈及目標在同一位 置，即r tm( )f r tt( )f 。如上所提，這個條件必頇經過準確的中途及終端導引才有 可能成立。因此探討攔截的可能性我們考慮其必要條件。其中一個最明顯的必要 條件是 ( ) ( ) m f t f r t  r t (4.1) 其中 . 代表向量的長度。如果(4.1)式不滿足，攔截一定無法達成；如果(4.1)式 滿足，則加上適當的導引律攔截有達成的可能。如圖 4.1 所示，當(4.1)式成立表 示當時間t_f時，飛彈與目標同時位於以圓心為座標原點的某個圓球上。

圖 4.1：目標與飛彈對座標原點的距離示意圖 除了(4.1)式的條件之外，成功攔截的必要條件還包括 ( ) tz f r t h   (4.2) 1 2 f c h t    t t t t (4.3) min c c t t (4.4) t_d 0 (4.5) ( ) mz f r t h   (4.6) ( ) ( ) tz mz r r t r t h    (4.7) 上列所述六個必要條件的意義分別說明如下： (4.2)式代表攔截的目標高度必頇高於 h 公里，此為確保截殺彈頭在空氣足夠稀薄 的高空運作簡化姿態控制；(4.3)式表示攔截時飛彈必頇是在終端導引階段，其中 h t 為第二節火箭脫離至目標鎖定所需的最短時間，為一個定值；(4.4)式表示滑行

時間必頇大於一個固定的最小值，以確保第一節火箭有足夠的時間可以安全脫 離；(4.5)式要求飛彈的發射必頇在雷達偵測到目標之後；(4.6)式代表攔截的飛彈 高度必頇高於 h 公里；(4.7)式代表目標與飛彈的高度差必頇小於h 公里。 _r 在此將考慮比較簡單的情況下即t 為一固定值，所以飛彈的位置為時間 t 的_c 函數而目標的位置為時間 t 及延遲時間t 的函數，分別如(3.10)及(3.2)式所列。令_d 變數x₁ t_d，x2 t，且定義函數 2 2 1 2 ( , ) _t( ) _m( ) h x x  r t  r t 1( ,1 2) tz( ) g x x r t h 2( ,1 2) 1 c 2 h 2 g x x     t t t t x 3( ,1 2) 1 g x x  x 4( ,1 2) mz( ) g x x r t h 5( ,1 2) tz( ) mz( ) r g x x  r t r t h 則(4.1)、(4.2)、(4.3)、(4.5)、(4.6)及(4.7)的必要條件可分別表示為 1 2 ( , ) 0 h x x  (4.8) 1( ,1 2) 0 g x x  (4.9) 2( ,1 2) 0 g x x  (4.10) 3( ,1 2) 0 g x x  (4.11) 4( ,1 2) 0 g x x  (4.12) 5( ,1 2) 0 g x x  (4.13) 其中(4.8)為等式限制式，(4.9)至(4.13)為不等式限制式。

4.2 最佳化問題的建立

令函數 1 2 1 2 ( , ) f x x  x x (4.14) 即 f x x( ,_{1 2}) t_d t，為雷達偵測到目標至飛彈發射後 t 秒的總時間。我們考慮的 最佳化問題之一為 (P1) 1 2 1 2 , max ( , ) x x f x x 1 2 1 1 2 3 1 2 4 1 2 ( , ) 0 ( , ) 0 ( , ) 0 ( , ) 0 h x x g x x subject to g x x g x x    _   _   _  透過此最佳化問題可以求出從雷達偵測到目標至攔截時間之最大值與其所 對應的延遲時間t ，飛彈飛行時間_d t_f，及預估的攔截點r t_t( )_f 。所求得的x 及1 x ，2 仍必頇代入檢查(4.10)限制式是否滿足，即g x x2( ,1 2)0，如果(4.10)式不滿足即 代表飛彈尚未進入終端階段不可能攔截成功，飛彈不應發射。這樣的狀況通常發 生於目標被偵測時已經太接近雷達站(或飛彈的發射點)。如果(4.10)式滿足則將 所求得的延遲時間t 視為容許最大的延遲時間。 _d

我們考慮的另一個最佳化問題是 (P2) 1 2 1 2 , min ( , ) x x f x x 1 2 1 1 2 2 1 2 3 1 2 4 1 2 5 1 2 ( , ) 0 ( , ) 0 ( , ) 0 ( , ) 0 ( , ) 0 ( , ) 0 h x x g x x g x x subject to g x x g x x g x x    _      _   _    透過此最佳化問題可以求出從雷達站偵測到目標至攔截時間之最小值與其 所對應的延遲時間t ，飛彈飛行時間_d t_f，及預估的攔截點r t_t( )_f 。與問題(P1)不同 的是，(P2)求取 f x x 的最小值且限制式包括(4.10)式。只有在(P1)有解的情況( ,_{1 2}) 下才考慮(P2)問題，並將所求得的延遲時間t 視為最短的延遲時間。 _d

4.3 延遲時間最大值與最小值

如果最佳化問題(P1)及(P2)均有可行解，且定義其分別算得的延遲時間為 max d

t 及tdmin，且0tdmin tdmax ，則時間區間



tdmin,tdmax



即為延遲時間的窗口，將

可提供攔截飛彈系統參考。大致上來說，雷達偵測到目標之後，飛彈至少需等待

min d

t 秒之後才能發射且不能等待超過t_dmax秒，以確保攔截有可能成功。我們可以

狀況 發射前的延遲時間t _d 攔截結果

Case 1 td tdmin 失敗

Case 2 tdmin  td tdmax 有可能成功

Case 3 td tdmax 失敗

第五章

最佳化演算法的發展

根據第四章所討論的最佳化問題可以用下列標準非線性規劃的問題表示： ( ) Minimize f x ( ) 0, 1, , ( ) 0, 1, , j i h j n subject to g i m        x x _(5.1) 其中 l x ，g_i( )x 0為不等式限制式，h_j( )x 0為等式限制式。一般來說，除 了非常特殊且簡單的情況，(5.1)的標準問題並沒有解析解，而必頇使用遞迴的方 式求解。因此在本章先後介紹兩種演算法來計算問題(5.1)的解，第一種方法是處 罰函數(Penalty Function，PF)的演算法。利用所謂的處罰函數是嘗詴將含有等式 及不等式的最佳化問題轉換成對等的沒有限制式(unconstrained)的問題，以方便 一般求最小值方法的使用，例如梯度法(gradient method)。第二種方法是所謂的 連續二次規劃(Sequential Quadratic Programming，SQP)的演算法。並設計一個簡 單的非線性最佳化問題來測詴所提出的兩種演算法是否可行，最後再比較兩種方 法的差異及可行性。

5.1 處罰函數(PF)演算法

考慮問題(5.1)定義處罰函數





2 2 1 1 ( ) max(0, ( )) ( ) m n i j i j P g h   







x x x (5.2) 及新的目標函數 ( , )c f( ) cP( )  x  x  x (5.3)

其中 c 為正數。PF 方法的主要精神是如果(5.1)的限制式不滿足，則(5.2)的P x( )大 於零，如果參數 c 很大，則將(5.3)式最小化可確定P x( )很小，當 c 時，則 ( ) 0 P x  ，即限制式均滿足。據此將(5.1)的問題轉換成下面沒有限制式的問題



( , ), l



Minimize  c x x (5.4)

在此利用 Gradient Search method，求出 x 使得( , )c x 最小。數值求解(5.4)需

要選擇三個參數分別為成長函數1，停止參數 0，及初始處罰參數c₀ 0。 定義梯度向量為 1 2 ( , ) ( , ) ( , ) ( , ) x l c x c x c c x         _        _   _{  }   _         x x x x 詳細求解步驟如下： ===================================================================== 步驟 1：選擇，，c ，及初始值₀ 0 x 步驟 2：迴圈求解，k 1, 2,3, 由 k1 x 開始尋找 k x 使得( ,c_k xk)最小，並記錄下每個 k x 所對應的( ,c_k xk) 步驟 3：計算 1 1 ( , k ) x ck    x 。 步驟 4： 1 1 1 ( , ) k k k k x ck         x x x 其中k 為搜尋距離，k 可利用 ( 1, ) 0 k k d c d    x 求得。

以下過程為判斷是否達到收斂條件： If ( xkxk1  or f(xk) f(xk1)  ) 跳出迴圈， k x 即為預估最小值的最佳解； Else c_k   c_k1，返回步驟 2，重新產生( ,c_k xk)； End =====================================================================

5.2 PF 範例

經過 5.1 節詳細介紹 PF 演算法後，接下來要透過 Matlab 程式撰寫求解一個 非線性的最佳化問題，並用圖形來驗證所求出的最小值與最大值是否接近實際的 最小值與最大值。 5.2.1 PF 範例最小值問題 1 2 ( ) Minimize f x  x x 2 2 1 2 2 1 2 ( ) 3 2 2 0 ( ) 2 3 2 0 h x x subject to g x x             x x 考慮本題之最小值問題，選擇1.5， 0.01，c0 0.05，初始值設定為 1 2 ( ,x x )  ( 1, 2)，經 Matlab 程式執行共 11 次後收斂，當( ,x x_{1 2}) ( 0.4865, 0.8200) ， 目標函數有最小值 f x( ) 1.3065。表 5.1 為 PF 求解結果的相關數據。 k x 1 x 2 c k f x x ( ,1 2) 0 -1.0000 -2.0000 0.05 -3.0000 1 -0.5873 -1.2205 0.05 -1.8078

2 -0.6492 -1.0172 0.075 -1.6664 3 -0.6426 -0.9662 0.1125 -1.6088 4 -0.6088 -0.9198 0.1688 -1.5286 5 -0.5776 -0.8845 0.2531 -1.4621 6 -0.5519 -0.8600 0.3797 -1.4119 7 -0.5312 -0.8436 0.5695 -1.3748 8 -0.5149 -0.8330 0.8543 -1.3479 9 -0.5024 -0.8264 1.2814 -1.3288 10 -0.4931 -0.8224 1.9222 -1.3155 11 -0.4865 -0.8200 2.8833 -1.3065 表 5.1：PF 求解範例最小值相關數據 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 x₁ x2 A (min) h (x)=0 B g (x)=0 圖 5.1：PF 範例最小值驗証 我們可以由圖 5.1 來檢視利用 PF 作法是否接近實際的最小值。其中綠色線 段代表限制式g( )x 0，綠色線段下方代表g( )x 0之解集合，藍色線段代表限

制式h( )x 0，因此滿足等式及不等式限制式之可行解區域為劣弧 AB。利用 PF 演算法所求出的可行解為 A 點，即( ,x x_{1 2}) ( 0.4865, 0.8200) 。紅色線段代表當 目標函數最小值為-1.3065，剛好通過 A 點，而黑色線段則代表目標函數值為 -1.5000 時，若求出的可行解在黑色線段上，則不同時滿足g( )x 0與h( )x 0。 因此利用 PF 演算法可以有效求得 A 點已經相當接近實際的最小值。 5.2.2 PF 範例最大值問題 1 2 ( ) Maximize f x  x x 2 2 1 2 2 1 2 ( ) 3 2 2 0 ( ) 2 3 2 0 h x x subject to g x x             x x 考慮本題之最大值問題，則選擇1.5， 0.01，c₀ 0.05，初始值設定為 1 2 ( ,x x )(1,1)，經 Matlab 程式執行共 11 次後收斂，當( ,x x1 2)(0.4913, 0.8042) ， 目標函數有最大值 f x( ) 0.3129。表 5.2 為 PF 求解結果的相關數據。 k x 1 x 2 c k f x x ( ,1 2) 0 1.0000 1.0000 0.05 2.0000 1 0.5025 0.6822 0.05 1.1847 2 0.3994 0.2719 0.075 0.6713 3 0.8037 -0.2417 0.1125 0.5620 4 0.5807 -0.4299 0.1688 0.1508 5 0.6466 -0.6633 0.2531 -0.0167 6 0.5708 -0.7105 0.3797 -0.1397 7 0.5540 -0.7562 0.5695 -0.2022

8 0.5259 -0.7744 0.8543 -0.2485 9 0.5123 -0.7897 1.2814 -0.2774 10 0.4987 -0.7977 1.9222 -0.2990 11 0.4913 -0.8042 2.8833 -0.3129 表 5.2：PF 求解範例最大值相關數據 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 x1 x2 h (x)=0 B (max) g (x)=0 A 圖 5.2：PF 範例最大值驗証 同樣地，我們可以由圖 5.2 來檢視利用 PF 作法是否接近實際的最大值。 其中綠色線段代表限制式g( )x 0，綠色線段下方代表g( )x 0之解集合，藍色 線段代表限制式h( )x 0，因此滿足等式及不等式限制式之可行解區域為劣弧 AB。利用 PF 演算法所求出的可行解為 B 點，即( ,x x_{1 2})(0.4913, 0.8042) 。紅 色線段代表當目標函數最大值為-0.3129，剛好通過 B 點，而黑色線段代表目標 函數值為 0.1 時，若求出的可行解在黑色線段上，則不同時滿足g( )x 0與 ( ) 0 h x  。因此利用 PF 演算法可以有效求得 B 點已經相當接近實際的最大值。

5.3 連續二次規劃(SQP)演算法

連續二次規劃方法 SQP，是常見且最簡單的非線性規劃最佳化的求解工具。 考慮標準問題(5.1)，以下說明僅考慮限制式為等式部份。 定義等式限制式為 1( ) ( ) ( ) n h h            x h x x 使用 Lagrange 函數為 1 ( , ) ( ) ( ) n j j j L f  h   



x λ x x (5.5) 其中_j為 Lagrange multipliers。 根據KKT(Karush-Kuhn-Tucker)必要條件為 1 ( ) ( ) n j j j f  h   x 



 x 0 (5.6) ( ) 0 j h x  (5.7) 即 ( , ) ( , ) xL L      _    x λ 0 x λ 。針對(5.6)式及(5.7)式左邊的函數利用泰勒式在 ( , ) k k x λ 展開可 得到 2 ( , ) ( , ) ( ) ( ) ( , ) k k k k k x x k T k k L L L            _   _           x λ x λ h x x 0 λ h x 0 x λ (5.8) 利用(5.8)式計算x 及λ 並令 1 k _ k _ x x x (5.9) 1 k  k  λ λ λ (5.10) (5.8)式、(5.9)式、(5.10)式所列的三個方程式為主要的遞迴步驟。

其中，在(5.8)式中 Lagrange 函數的 Hessian 矩陣為 2 2 2 1 ( , ) ( ) ( ) n x j j j L f  h   x λ   x 



 x (5.11) 將(5.8)式經過移項整理之後可以得到下列矩陣 2 _{( ) ( )} ( , ) ( ) ( ) ( ) k k k k k k k T k f L                _    _       x x λ h x x x h x λ λ h x 0 h x (5.12) 在此假設 k  d x，並將(5.10)式代入(5.12)式後可得到下列矩陣 2 1 ( ) ( ) ( , ) ( ) ( ) ( ) k k k k k k k k T k k k f L                  _          x x λ h x d x h x λ h x 0 λ λ h x 2 1 ( ) ( ) ( ) ( , ) ( ) ( ) ( ) k k k k k k k k k k T k k f L                     _         x x λ h x d x h x λ h x λ h x 0 λ h x 2 1 ( ) ( , ) ( ) ( ) ( ) k k k k k x k T k k f L                _        d x x λ h x h x 0 λ h x (5.13) 上述(5.13)式稱為 Newton-Lagrange System (LNS)。 定義不等式限制式為 1( ) ( ) ( ) m g g            x g x x 在此考慮不等式限制式gi( )x 0與等式限制式hj( )x 0同時存在時，則只有 當不等式限制式gi( )x 0時，才將此不等式限制式加入 LNS 中，其矩陣為 2 1 1 ( , , ) ( ) ( ) ( ) ( ) 0 0 ( ) ( ) 0 0 ( ) k k k k k k k x k T k k k T k k L f            _   _{ }         _            x λ β h x g x d x h x λ h x g x β g x (5.14)

其中 Lagrange 函數改變為 1 1 ( , , ) ( ) ( ) ( ) n m j j i i j i L f  h g    







x λ β x x x (5.15) 定義收斂參數： 0 詳細求解步驟如下： ===================================================================== 步驟 1：決定初始值 0 0 0 (x λ β, , ) 步驟 2：迴圈求解 LNS，k1, 2,3, ，可求出新的增加量 k d 、_{λ 與}k1 k1 β 步驟 3：更新 k1  k  k x x d 步驟 4：檢查是否收斂，收斂條件 kT  k  d d ，若滿足，則停止。 步驟 5：令k  k 1，重新回到步驟 2 直至收斂為止 =====================================================================

5.4 SQP 範例

經過 5.3 節詳細介紹 SQP 演算法後，接下來要透過 Matlab 程式撰寫求解一 個非線性的最佳化問題，此範例與 5.2 節所提出的 PF 範例相同，目的是要來比 較所求出的可行解是否接近，最後使用表 5.5 及表 5.6 來比較兩種演算法的差異 及可行性。 5.4.1 SQP 範例最小值問題 1 2 ( ) Minimize f x  x x 2 2 1 2 2 1 2 ( ) 3 2 2 0 ( ) 2 3 2 0 h x x subject to g x x             x x

若考慮本題之最小值問題，選擇 0 0 ( , )(0.5, 0.5)， 6 10  _  ，初始值設定 為( ,x x_{1 2})(2, 5) ，經 Matlab 程式執行共 7 次後收斂，當 1 2 ( ,x x ) ( 0.4725, 0.8155) ，目標函數有最小值 ( )f x  1.2880。 k x 1 x 2 dk 1 k   k1 1 2 ( , ) f x x 1 0.6022 -2.8387 (-1.3978,2.1613) 0.2661 0.5000 -2.2365 2 -0.4373 -1.8304 (-1.0394,1.0083) 0.1826 0.5000 -2.2677 3 -0.8183 -0.9735 (-0.3810,0.8569) 0.2221 0.5000 -1.7918 4 -0.5510 -0.8215 (0.2672,0.1520) 0.3417 0.0651 -1.3725 5 -0.4781 -0.8155 (0.0729,0.0059) 0.3337 0.0295 -1.2936 6 -0.4726 -0.8155 (0.0056,9.0850e-06) 0.3333 0.0291 -1.2881 7 -0.4725 -0.8155 (3.2821e-05,2.1267e-11) 0.3333 0.0291 -1.2880 表 5.3：SQP 求解範例最小值相關數據 5.4.2 SQP 範例最大值問題 1 2 ( ) Maximize f x  x x 2 2 1 2 2 1 2 ( ) 3 2 2 0 ( ) 2 3 2 0 h x x subject to g x x             x x 若考慮本題之最大值問題，則選擇 0 0 ( , )(0.5, 0.5)， 6 10  _  _{，初始值設定為} 1 2 ( ,x x )(2, 5) ，經 Matlab 程式執行共 7 次後收斂，當( ,x x_{1 2})(0.4725, 0.8155) ， 目標函數有最大值 ( )f x  0.3430。

k x 1 x 2 dk 1 k   k1 1 2 ( , ) f x x 1 1.4624 -2.3226 (-0.5376,2.6774) 0.2177 0.5000 -0.8602 2 1.7495 -0.4148 (0.2872,1.9078) 0.0712 0.5000 1.3347 3 0.9609 -0.8677 (-0.7886,-0.4529) 0.0195 0.3871 0.0932 4 0.5972 -0.8162 (-0.3637,0.0515) 0.0323 0.3694 -0.2190 5 0.4855 -0.8155 (-0.1116,6.8202e-04) 0.0631 0.4020 -0.3300 6 0.4727 -0.8155 (-0.0128,1.1985e-07) 0.0752 0.4152 -0.3428 7 0.4725 -0.8155 (-1.7424e-04,3.7251e-15) 0.0757 0.4156 -0.3430 表 5.4：SQP 求解範例最大值相關數據 由於本題與 5.2 節 PF 範例完全相同，因此為了檢視 SQP 與 PF 演算法的差 異及可行性，可參考表 5.5 及表 5.6 分別為利用兩種方法來求解最小值及最大值 的求解結果。 方法 執行次數 可行解x ₁ 可行解x ₂ 目標函數最小值 PF 11 -0.4865 -0.8200 -1.3065 SQP 7 -0.4725 -0.8155 -1.2880 表 5.5：最小值問題方法比較 方法 執行次數 可行解x ₁ 可行解x 2 目標函數最大值 PF 11 0.4913 -0.8042 -0.3129 SQP 7 0.4725 -0.8155 -0.3430 表 5.6：最大值問題方法比較

透過表 5.5 及表 5.6 中可以發現兩種方法所求得的可行解相當接近，且由圖 5.1 及圖 5.2 可得知求出的可行解均相當接近實際的最小值與最大值。經 Matlab 程式執行以 SQP 方法執行次數較少，PF 執行次數較多但是最終亦可準確求出可 行解。因此若不考慮執行次數多寡，PF 與 SQP 演算法皆可作為我們所要探討最 佳化問題(P1)及(P2)的求解工具。

第六章

數值計算及模擬

6.1 飛彈的推力及物性

本論文飛彈的推力設定為修改探空五號[9]的資料。圖 6.1 分別為攔截飛彈在 飛行過程的質量變化與推力大小的設定。表 6.1 為比較各節火箭推力與作用時 間，表 6.2 為比較各節的燃料質量。表 6.3 是飛彈各節質量的分配狀況。 0 50 100 150 200 250 300 0 500 1000 1500 2000 Time (sec) M a s s ( k g ) 0 50 100 150 200 250 300 0 5 10 15x 10 4 Time (sec) T h ru s t (N ) 圖 6.1：飛行過程中的質量與推力隨時間變化示意圖 本論文火箭 探空五號 推進器 第一節 第二節 第一節 第二節 總衝量 161,500(kg-s) 854,20(kg-s) 102,968(kg-s) 148,560(kg-s) 總燃燒時間 16s 36s 5.4s 29s 表 6.1：比較探空五號與本論文火箭的推力與作用時間

本論文火箭 探空五號 第一節 703 (kg) 429.20(kg) 第二節 373 (kg) 645.17(kg) 燃料總質量 1076 (kg) 1074.37(kg) 表 6.2：比較探空五號與本論文火箭的燃料質量 飛彈總質量 1650(kg) 飛彈燃料總質量 1076(kg) 第一節燃料質量 703(kg) 第二節燃料質量 373(kg) 第一節燃料用完時質量 947(kg) 第一節飛彈脫離時質量 810(kg) 第二節燃料用完時質量 437(kg) 第二節飛彈脫離時質量 21(kg) 表 6.3：飛彈燃料質量分配狀況

6.2 模擬初始條件

分別建立三種不同來襲目標的資料庫如下表所示 初始位置(km) 初始速度(km/s) 最大高度(km) 目標一 (900,-900,-30) (-2.1,2.0,-2.0) 239.13 目標二 (750,-710,-30) (-1.8,1.9,-1.6) 161.65 目標三 (350,-280,-97) (-1.8,1.4,-0.1) 97.526 表 6.4：目標資料庫

目標一在雷達站的北方 900 公里，西方 900 公里，高度 30 公里，距離雷達 站為 1273.1 公里。目標二在雷達站的北方 750 公里，西方 710 公里，高度 30 公 里，距離雷達站為 1033.2 公里。目標三在雷達站的北方 350 公里，西方 280 公 里，高度 97 公里，距離雷達站為 458.59 公里。三個來襲目標均朝向東南方飛行 至雷達站附近落地，如圖 6.2 所示。 -1 0 1 2 3 4 5 6 7 8 9 x 105 -10 -8 -6 -4 -2 0 x 105 0 0.5 1 1.5 2 2.5 x 105 North (m) Target trajectory East (m) A lt it u d e ( m ) Target1 Target2 Target3 r_t3(-t_d) r_t1(-t_d) r_t2(-t_d) 圖 6.2：不同目標的飛行軌跡圖(NED) 圖 6.3 為三個來襲目標的飛行高度隨時間變化關係圖，可以很明顯看出目標 一飛行高度最高且被雷達偵測時處於逐漸上升階段，而目標三飛行高度最低且被 雷達偵測時處於正在下降階段，目標二飛行最大高度則介於目標一與目標三之間 且被雷達偵測時屬於逐漸上升階段。

0 50 100 150 200 250 300 350 400 450 0 0.5 1 1.5 2 2.5x 10 5 Time (sec) T a rg e t a lt it u d e ( m ) Target1 Target2 Target3 圖 6.3：不同目標的飛行高度對時間圖 飛彈初始位置(-1,0,-0.01) km。發射點在雷達站的南方 1 公里，高度 10 公尺。 飛彈靜止於發射點，發射仰角固定為 70 度，發射的方向角固定為粗估的攔截點 方向，即目標初始位置，更精確的發射方向角將在 6.4 節介紹。 其他時間參數t1 16s，tc 18s，t2 36s，th 10s。第一節飛彈脫離後的 質量m₂₀810kg，第二節燃燒時的最大推力T₂ 23580N，第二節推力的質量燃 燒速率 ₂ 373 36 m  kg/s，重力加速度g (0, 0, 0.0098)km/s2。攔截必要條件(4.9)式及 (4.12)式參數h30，(4.13)式參數hr 10。

根據 3.3 節所推導的目標位置依三個目標分別為 目標一 2 900 2.1 0 1 ( ) 900 2.0 ( ) 0 ( ) 2 30 2.0 0.0098 t d d r t t t t t               _{  } _   _{ }              ，t t_d 目標二 2 750 1.8 0 1 ( ) 710 1.9 ( ) 0 ( ) 2 30 1.6 0.0098 t d d r t t t t t               _{  } _   _{ }              ，t t_d 目標三 2 350 1.8 0 1 ( ) 280 1.4 ( ) 0 ( ) 2 97 0.1 0.0098 t d d r t t t t t               _{  } _   _{ }              ，t t_d 目標速度依三個目標分別為 目標一 2.1 0 ( ) 2.0 0 ( ) 2.0 0.0098 t d v t t t          _{  } _          目標二 1.8 0 ( ) 1.9 0 ( ) 1.6 0.0098 t d v t t t          _{  } _         

目標三 1.8 0 ( ) 1.4 0 ( ) 0.1 0.0098 t d v t t t          _{  } _          透過 Simulink 程式模擬將飛彈發射自由飛行至第二節點火前，根據 3.4 節所推導 飛彈的位置如下。當第一節火箭燃畢，即tt₁時， 飛彈的位置為 1 1 1.7311 ( ) 2.6115 6.2311 m r t R        _{ }     飛彈的速度為 1 1 0.3781 ( ) 0.3616 0.7866 m v t V        _{ }     當滑行完畢後，即t t₁ t_c時， 2 1 1.7311 0.3781 0 1 ( ) 2.6115 0.3616 18 0 18 2 6.2311 0.7866 0.0098 m c r t t               _{  }  _  _{ }             當第二節火箭點火之後，即t1    tc t t1 tc t2時，推力產生的加速度為 4 4 1 4 0.0133 1.6963 10 ( ) 0.0127 ( ) 1.6221 10 0.0214 2.7372 10 m c a t t t t              _{ }   _  _           飛彈的位置為 4 2 2 3 4 4 1.7311 0.3781 0 0.0133 1.6963 10 1 1 1 ( ) 2.6115 0.3616 ( 16) 0 ( 16) ( 34) 0.0127 ( 34) 1.6221 10 2 2 6 6.2311 0.7866 0.0098 0.0214 2.7372 10 m r t t t t t                          _{  }  _   _{ }    _ _  _  _                      

第二節火箭燃畢時，將t  t₁ t_c t₂代入上式，可得飛彈的位置為 4 2 2 3 4 2 4 1.7311 0.3781 0 0.0133 1.6963 10 1 1 1 2.6115 0.3616 54 0 54 36 0.0127 36 1.6221 10 2 2 6 6.2311 0.7866 0.0098 0.0214 2.7372 10 32.3895 m f r                          _{  }  _  _{ }    _{ }   _  _                        31.9299 50.9403           飛彈的速度為 4 2 4 2 4 0.3781 0 0.0133 1.6963 10 1 0.3616 0 54 36 0.0127 36 1.6221 10 2 0.7866 0.0098 0.0214 2.7372 10 0.9847 0.9417 1.2362 m f v                      _{  } _   _{ }   _  _                        _{ }     當t  t₁ t_c t₂時，即第三節截殺彈頭飛行時，飛彈的位置為 2 2 2 2 0 1 ( ) ( 70) 0 ( 70) 2 0.0098 36.5412 0.9847 33.9880 0.9417 59.6039 1.9222 0.0049 m m f m f r t r v t t t t t t          _{ }            _  _      

6.3 模擬結果

先透過 Simulink 程式實際模擬有中途導引狀況下當三種不同來襲目標與飛 彈交錯時的攔截狀況，之後再利用 PF 及 SQP 演算法求出延遲時間最大值及最小 值來驗證是否接近實際情況。

6.3.1 Simulink 模擬目標一來襲 利用 6.2 節所設定的初始條件，並且透過 Simulink 程式來模擬目標一在不同 延遲時間的攔截狀況如表 6.5 所示，可以得到當目標與飛彈交錯時的零力誤失距 離、目標高度、飛彈飛行時間、 ( )r t 與 ( )_m r t ，其中延遲時間每隔 10 秒模擬一次。 t 延遲時間 （s） 零力誤失 距離（km） 目標高度 (km) 飛行時間 （s） ( ) m r t (km) ( ) t r t (km) 0 123.48 182.94 323.93 343.00 387.90 10 102.88 174.73 321.70 327.64 364.37 20 82.266 166.08 319.43 311.80 340.91 30 67.998 160.40 314.28 302.35 326.16 40 54.313 154.76 308.88 293.20 311.94 50 41.479 149.31 303.18 284.50 298.58 60 29.470 144.08 297.15 276.35 286.08 70 18.688 139.25 290.84 269.82 274.78 80 15.745 139.11 280.78 269.47 274.47 90 13.374 138.84 271.06 269.87 273.84 100 11.276 138.02 261.66 268.74 271.95 110 9.4430 136.78 252.54 266.63 269.11 120 7.9499 135.13 243.70 263.54 265.33 130 6.6508 133.31 235.05 259.67 261.22 140 5.5609 131.20 226.09 255.09 256.47 150 4.6438 128.64 218.31 249.87 250.76 160 3.8183 125.91 210.22 244.15 244.71 170 3.0696 122.99 202.26 237.82 238.34

180 2.5609 119.74 194.42 231.09 231.31 190 2.0304 116.29 186.73 223.80 223.93 200 1.6978 112.64 179.11 216.27 216.22 210 1.3078 108.78 171.63 208.21 208.17 220 1.0306 104.70 164.25 199.77 199.76 230 0.9190 100.39 156.94 191.11 191.01 240 0.7070 95.858 149.72 182.12 181.91 250 0.6000 91.247 142.59 172.69 172.78 260 0.5168 86.223 135.55 163.11 162.98 270 0.4759 81.114 128.59 153.12 153.16 280 0.4123 75.569 121.70 143.04 142.66 290 0.3141 69.926 114.90 132.56 132.15 300 0.3090 64.185 108.14 121.84 121.63 310 0.4850 58.162 101.49 110.76 110.81 320 0.4321 51.661 94.980 99.374 99.363 330 0.8890 44.853 88.546 87.776 87.685 340 1.3120 37.346 82.403 75.425 75.233 表 6.5：模擬目標一在不同延遲時間的攔截情形 圖 6.4 為目標一在延遲時間為零秒的情況下，目標與飛彈的飛行軌跡圖，可 以看出當兩者交錯時飛彈的高度已明顯低於目標高度，造成最後飛彈比目標提早 落地而導致攔截失敗。目標一與飛彈詳細的飛行高度可參考圖 6.5。

-5 0 5 10 x 105 -10 -5 0 x 105 0 0.5 1 1.5 2 2.5 x 105 North (m) Missile & Target trajectory

East (m) A lt it u d e ( m ) Missile Target 圖 6.4：飛彈與目標一的飛行軌跡圖(NED) 圖 6.5 為目標一在延遲時間為零秒的情況下，飛彈飛行 367.60 秒之後比目標 提早落地，此時目標仍然在空中飛行且高度為 129.79 公里，導致最終攔截失敗， 因此當目標一來襲時必頇要有延遲時間等待目標接近才允許發射飛彈。 0 50 100 150 200 250 300 350 400 450 0 0.5 1 1.5 2 2.5x 10 5 Time (sec) A lt it u d e ( m ) Missile Target 圖 6.5：飛彈與目標一的飛行高度對時間圖

圖 6.6 為延遲時間為 70 秒至 340 秒的零力誤失距離圖，當延遲 70 秒發射時， 飛彈與目標一交錯的零力誤失距離為 18.69 公里，隨著延遲時間的增加，所對應 的零力誤失距離愈來愈小。 圖 6.7 同樣為延遲時間為 70 秒至 340 秒的目標高度對時間圖，同樣隨著延 遲時間的增加，所對應的目標高度愈來愈低，當延遲時間為 340 秒的目標高度接 近 30 公里，且由表 6.5 查表可知此時飛彈飛行總時間為 82.40 秒，剛進入第三節 終端飛行階段。 100 150 200 250 300 350 0 2 4 6 8 10 12 14 16 18 20 Delay Time (s) Z E M ( k m ) 圖 6.6：目標一不同延遲時間對應的零力誤失距離圖

100 150 200 250 300 20 40 60 80 100 120 140 Delay Time (s) T a rg e t A lt it u d e ( k m ) 圖 6.7：目標一不同延遲時間對應的高度關係圖 6.3.1.1 執行 PF 演算法求解目標一 當目標一來襲時，利用 5.1 節所發展的 PF 演算法，透過 Matlab 程式撰寫可 求出第四章所提出的最佳化問題(P1)及(P2)。圖 6.8 為根據第四章所提到的等式 限制式h x x( ,1 2)0及不等式限制式g x x1( ,1 2)0與g x x5( ,1 2)0所畫的示意圖，橫 軸代表延遲時間x₁0，縱軸代表飛彈飛行時間x₂ 0，其中藍色曲線代表 1 2 ( , ) 0 h x x  ；綠色曲線代表g x x1( ,1 2)0，綠色曲線下方即為g x x1( ,1 2)0的可行 解區域；紅色曲線代表g x x₅( ,_{1 2})0，紅色曲線上方即為g x x₅( ,_{1 2})0的可行解區 域。因此從圖形可看出欲同時滿足等式及不等式限制式的延遲時間窗口為藍色線 段m m ，即延遲時間最小值為_{1 2} m₁ 115，延遲時間最大值為m₂ 321。

0 50 100 150 200 250 300 350 400 0 50 100 150 200 250 300 350 400 x₁ x2 X: 321.3 Y: 86.98 X: 115.8 Y: 232.9 g1(x1,x2)=0 h(x1,x2)=0 g₅(x₁,x₂)=0 m₂ m₁ 圖 6.8：目標一延遲時間示意圖 (i).欲求解max : ( ,f x x1 2) x1 x2，即延遲時間與飛行時間最大值 兩個初始值分別為( ,x x_{1 2})(335,90)及( ,x x_{1 2})(280,120)，詳細求解結果如表 6.6 所示 演算法 初始值 1 2 ( ,x x ) 延遲時間 最大值（s） 飛行時間 （s） 目標函數值 （s） PF (335,90) 320.95 87.189 408.13 (280,120) 320.39 87.550 407.94 表 6.6：PF 求解目標一延遲時間最大值