Airplane Model

We consider a networked control system to perform the turning maneuver of an airplane in which a main controller directs subcontrollers for the ailerons and the rudder.

To make a turn, an aircraft rolls towards the direction of the turn by moving its ailerons (surfaces attached to the end of the wings). The rolling causes a yawing moment in the opposite direction, called adverse yaw, which is countered by using its rudder (a surface attached to the vertical stabilizer).

The system consists of three controllers. The main controller determines the desired control angles, and the subcontroller (for the ailerons and the rudder) gradually change the surface angles according to the main controller.

Dynamics

We consider both nonlinear and linear systems of differential equations for this model. The lateral dynamics of an aircraft can be specified as nonlinear ordinary differential equations, where \(\beta\) is the yaw angle, \(p\) is the rolling moment, \(r\) is the yawing moment, \(\phi\) is the roll angle, \(\psi\) is the direction angle, and \(Y(\beta,\delta_a,\delta_r)\), \(L(\beta,\delta_a,\delta_r)\), and \(N(\beta,\delta_a,\delta_r)\) are (linear) functions denoting the coefficients of side force, rolling moment, and yawing moment.

\begin{aligned} \frac{d\beta}{dt} = & \frac{Y(\beta,\delta_a,\delta_r)}{m V} - r + \frac{V}{g} \cos(\beta) \sin(\phi) \\
\frac{dp}{dt} = & (c_1 r + c_2 p)\cdot r \tan(\phi) + c_3 L(\beta,\delta_a,\delta_r) + c_4 N(\beta,\delta_a,\delta_r) \\
\frac{dr}{dt} = & (c_8 p - c_2 r) \cdot r \tan(\phi) + c_4 L(\beta,\delta_a,\delta_r) + c_9 N(\beta,\delta_a,\delta_r) \\
\frac{d\phi}{dt} = & p \\
\frac{d\psi}{dt} = & \frac{g}{V} \tan(\phi) \\
\end{aligned}

These nonlinear ODEs can be approximated as the linear ODEs, assuming small perturbations:

\begin{aligned} \frac{d\beta}{dt} = & Y_\beta \beta - r + \frac{V}{g} \phi + Y_r \delta_r \\
\frac{dp}{dt} = & L_b \beta + L_p p + L_r r + L_a \delta_a + L_r \delta_r \\
\frac{dr}{dt} = & N_b \beta + N_p p + N_r r + N_a \delta_a + N_r \delta_r \\
\frac{d\phi}{dt} = & p \\
\frac{d\psi}{dt} = & \frac{g}{V} \phi \\
\end{aligned}

Subcontrollers

In each round, a subcontroller \(M\) receives a goal angle \(\mathit{goal}_m\) from the main controller, determines the new moving rate \(\mathit{rate}_M\) based on \(\mathit{goal}_m\) and the current sampled value \(v_M\) of the surface angle \(\alpha_M\), and sends back the current angle \(v_M\) to the main controller. The new moving rate \(\mathit{rate}_M\) can be given by \( \mathrm{sign}(\mathit{goal}_M - v_M) \cdot \min(\mathrm{abs}(\mathit{goal}_M - v_M) / T, \mathit{max}_M)\), and the continuous dynamics of the surface angle is given by the ODE \(\frac{\mathrm{d}\alpha_M}{\mathrm{d}t} = \mathit{rate}_M\).

Main Controller

The main controller determines the goal angles for the subcontrollers to make a coordinated turn. In each round, it receives a desired direction \(\mathit{goal}_\psi\) of the airplane (from the pilot) and the surface angles \((v_L, v_V, v_R)\) from the subcontrollers, and then sends back the new goal angles \((\mathit{goal}_L,\mathit{goal}_V,\mathit{goal}_R)\) to the subcontrollers, based on the current sampled position values \((v_\psi, v_\phi, v_\beta)\) of the direction angle \(\psi\), the roll angle \(\phi\), and the yaw angle \(\beta\).

We consider a simple control logic to decide the new goal angles \((\mathit{goal}_L,\mathit{goal}_V,\mathit{goal}_R)\), by three control functions: \(f_\phi(v_\phi, v_\psi, \mathit{goal}_\psi)\) for the goal roll angle \(\mathit{goal}_\phi\), \(f_{R}(v_\phi, \mathit{goal}_\phi)\) for the goal angle \(\mathit{goal}_R\) of the right aileron (where \(\mathit{goal}_L\) for the left aileron is \(- \mathit{goal}_R\)), and \(f_V(v_\beta)\) for the goal rudder angle \(\mathit{goal}_V\). For example:

\begin{aligned} f_\phi(v_\phi, v_\psi, g_\psi) &= v_\phi + \mathrm{sign}(h - v_\phi) \cdot \min(\mathrm{abs}(h - v_\phi),1.5), &&\mbox{where}\;\; h = 0.32 (\mathit{g}_\psi - v_\psi), \\
f_{R}(v_\phi, g_\phi) &= \mathrm{sign}(\mathit{g}_\phi - v_\phi)\cdot \min(0.3 \,\mathrm{abs}(\mathit{g}_\phi - v_\phi), 15)), \\
f_V(v_\beta) &= \mathrm{sign}(-v_\beta)\cdot \min(0.3 \,\mathrm{abs}(v_\beta), 10). \end{aligned}

Bounded Reachability

We have first performed bounded reacability analysis for the model with simplified controls. We consider three variants: a single-rate nonlinear model where every controller has a 0.5 period, a single-rate approximated model by linear differential equations, and a multirate linear model where the aileron controller has a 0.25 period and the rudder controller has a 0.5/3 period. We set a timeout of 30 hours for the experiments.

Files

To generate SMT files for this model in the new/standard encodings, we have developed a simple python script. For example, we can run the following commands to generate the SMT file using the new encoding and check its satisfiability using dReal for the unsat case of the simplified airplane model with bound \(k = 2\), provided that the library file gen.py is in the same directory:

python airplane-single-p.py 2 > airplane-single-p_2.smt2
dReal airplane-single-p_2.smt2

For the dReach scripts, we use the option -k N -l N to set both lower and upper bound to \(N\). For example, for the simplified airplane dReach model, we can set \(N = 2\) using the following command:

dReach -k 2 -l 2 airplane-single.drh

The following are the python scripts (to generate SMT files) and the dReach models for the three airplane models.

• SMT generation script: gen.py
• Nonlinear model
  • Unsat: New encoding, Standard encoding, dReach script
  • Sat: New encoding, Standard encoding, dReach script
• Simplified model
  • Unsat: New encoding, Standard encoding, dReach script
  • Sat: New encoding, Standard encoding, dReach script
• Multirate model
  • Unsat: New encoding, Standard encoding, dReach scripts(k=4n, k=4n+2)
  • Sat: New encoding, Standard encoding, dReach scripts(k=4n, k=4n+2)

Compositional Bounded Reachability

We have applied compositional reasoning to conduct a bounded reachability analysis in a compositional way. Assuming that the goal angles chosen by the main controller are close to the current surface angles of the ailerons and the rudder (i.e., the differences are less than \(0.05\)), we have performed bounded model checking up to \(k = 20\) from the initial condition \((\beta, \phi, p, r, \delta_a, \delta_r) = \vec{0}\) for the safety property \(\lvert\beta\rvert \leq 0.2 \) only using the main controller, which took less than 2 minutes.

The assumptions are then proved by using compositional inductive reasoning. For example, for the left aileron, we show that if the difference between the goal \(g_L\) angle and the current angle \(v_L\) is less than \(0.03\) at the beginning of the period, then the difference between the goal angle \(g_L\) and the surface angle \(x_L\) is: (i) always less than \(0.05\) during the period, and (ii) less than \(0.03\) again at the end of the period, taking into account the maximal local clock skew \(\epsilon < 0.0001\), sampling times \(0.001\), and actuator response times \(0.003\) (which took about 2 hours).

Files

The following are the dReach scripts for the compositional bounded reachability analysis. The following commands are used for the analysis:

dReach -k 20 airplane-main.drh
dReach -k 2 airplane-ail.drh
dReach -k 2 airplane-ail-ind.drh

The first file contains the negation of the formula stating that \(\lvert\beta\rvert \leq 0.2 \) for the main controller, provided the assumptions. The other files are regarding the compositional inductive analysis of the left aileron (the other cases are similar). The second file contains the negation of the formula statings that if \( \lvert g_L - v_L \rvert < 0.03\) at the beginning of the period, then \( \lvert g_L - x_L \rvert < 0.05\) during the period. The third file contains the negation of the formula stating that if \( \lvert g_L - v_L \rvert < 0.03\) at the beginning of the period, then \( \lvert g_L - x_L \rvert < 0.03\) again at the end of the period.

• main compositional reachability
• left compositional inductive 1 (k = 2)
• left compositional inductive 2 (k = 2)