Missile Optimization Guidance

Problem Description

Generally, assuming the target is stationary or moving at low speed, the geometric representation of guidance within a two-dimensional horizontal plane is shown as follow:

Among them $V_i$ are the velocities of the $i$ th missile, respectively. The remaining flight time of the i-th missile is represented by guidance flight state parameters such as line of sight distance $r_i$ , line of sight angle $q_i$, rake angle $\varphi_i$ , and trajectory angle $\theta_i$ , as shown in the figure (for simplicity, the index of the i-th missile is omitted) Missile guidance follows the so-called proportional guidance relationship, which is $\varphi=\theta − q, \dot{\theta} = N(t) \dot q$. Among them $N$ is the guidance gain (proportional guidance coefficient) that varies over time, Therefore, the dynamic equation of missile line of sight angle can be expressed as:

$$\begin{split} \dot{r}&= -V \cos \varphi \\ r \dot{q}&= -V \sin \varphi \\ r \dot{\varphi}&= (1-N(t)) V \sin \varphi \\ \dot{V}&= -p_0 V^2 - p_1 \cdot \tfrac{V^2 \sin^2 \varphi}{r^2} \cdot N(t) \end{split} \tag{1}$$

Among them $p_0, p_1$ are coefficients related to aerodynamic lift, aerodynamic drag, and missile parameters.

The lateral acceleration of the missile is described as:

$$a = V \dot{\theta} = V N \dot{q} = - \frac{N V^2 \sin \varphi}{r} \tag{2}$$

Simply put, consider the limitation of maximum lateral overload during missile guidance:

$$r a_{min} \leq - N V^2 \sin \varphi \leq r a_{max} \tag{3}$$

In addition, based on the constraints of landing angle and finite field of view, the guided flight process also includes the following boundary constraints:

$$\begin{split} &r(t_0) = r_0, \ q(t_0) = q_0, \ \varphi(t_0) = \varphi_0 \\ &\varphi_{min} \leq \varphi \leq \varphi_{max} \end{split} \tag{4}$$

The optimization objective comprehensively considers overload (least squares) integration and maximizing end velocity, and the OCP model is defined as:

$$\begin{array}{cll} \underset{r(\cdot), q(\cdot), \varphi(\cdot)}{\operatorname{minimize}} & -V(t_f) + \int_{0}^{T} \ell(a(t)) \mathrm{d} t \\ \text { subject to } & r(0)=\bar{r}_{0}, q(0)=\bar{q}_{0}, \varphi(0)=\bar{\varphi}_{0} & \\ & (1),(3),(4). \end{array}$$

Modeling

N = 100
prob = multi_stage_problem('guidance', NN + 1)
g = 9.81
t0 = 0
vMin, vMax = 1.0, 1000.0
ngMin, ngMax = -16.0, 16.0

r0, q0, v0, varphi0 = prob.parameters(['r0', 'q0', 'v0', 'varphi0'],  stage_dependent=False)
rf, qf, varphif = prob.parameters(['rf', 'qf', 'varphif'],  stage_dependent=False)

rMin, rMax = prob.parameters(['rMin', 'rMax'],  stage_dependent=False)
qMin, qMax = prob.parameters(['qMin', 'qMax'],  stage_dependent=False)
varphiMin, varphiMax = prob.parameters(['varphiMin', 'varphiMax'],  stage_dependent=False)
aMin, aMax = prob.parameters(['aMin', 'aMax'],  stage_dependent=False)
tfMin, tfMax = prob.parameters(['tfMin', 'tfMax'],  stage_dependent=False)

vfWeight, aWeight, ngWeight = prob.parameters(['vfWeight', 'aWeight', 'ngWeight'],  stage_dependent=False)

r = prob.variable('r', hard_lowerbound=rMin, hard_upperbound=rMax)
q = prob.variable('q', hard_lowerbound=qMin, hard_upperbound=qMax)
v = prob.variable('v', hard_lowerbound=vMin, hard_upperbound=vMax)
varphi = prob.variable('varphi', hard_lowerbound=varphiMin, hard_upperbound=varphiMax)
ng = prob.variable('ng', hard_lowerbound=ngMin, hard_upperbound=ngMax)

tf = prob.variable('tf', hard_lowerbound=tfMin, hard_upperbound=tfMax)

rho = 1.211
CD0 = 0.005
Sref = 0.29
m = 200.0
K = 0.03
p0 = 0.5 * rho * CD0 * Sref / m
p1 = 2 * m * K / (rho * Sref ** 2)

""" ode """
dr = -v * cos(varphi)
dq = -v * sin(varphi) / (r)
dv = -p0 * v ** 2 - p1 * (v * sin(varphi) / (r)) ** 2 * ng ** 2
dvarphi = (1 - ng) * v * sin(varphi) / (r)
dxdt = [dr, dq, dv, dvarphi, 0]
ts = tf / NN
ode = differential_equation(
    state=[r, q, v, varphi, tf],
    state_dot=dxdt, 
    stepsize=ts, 
    discretization_method='trapezoid')
prob.equality(ode)

# inequality
a = ng * v ** 2 * sin(varphi) / (r)
ineqa = general_inequality(expr=[a, a], sign=['>=', '<='], bound=[aMin, aMax])
prob.inequality(ineqa)

# minimize fuel consumption
obj_a = general_objective(aWeight * a ** 2 + ngWeight * ng ** 2)
prob.objective(obj_a)
obj_v = general_objective(-vfWeight * v)
prob.end_objective(obj_v)

seq = general_equality([r - r0, q - q0, v - v0, varphi - varphi0])
prob.start_equality(seq)

eeq = general_equality([r - rf, q - qf, varphi - varphif])
prob.end_equality(eeq)

option = codegen_option()
codegen = code_generator()
codegen.codegen(prob, option)

Solution

---

[1]. Jiang, Huan, et al. "Cooperative guidance with multiple constraints using convex optimization." Aerospace science and technology 79 (2018)