F-16 FCS Issues

[seanmcleod]

@dpculp commented in a response to another discussion #809 (comment) that he was seeing an issue with the JSBSim repo's F-16 FCS for a simple test case.

I just ran a cruise test with the F-16 model and it suffered a pitch excursion at about 90 seconds into cruise flight leading to a crash.

In this test the airplane is trimmed at 20k feet altitude and let loose for 5 minutes with no control inputs.

There have been reports in the past about pitch excursions with this model, and it's probably coming from integrator wind-up in a PID controller, although it could be coming from elsewhere in the pitch channel.

[seanmcleod]

Taking a look to see if I can reproduce the results of @dpculp.

import jsbsim
import matplotlib.pyplot as plt

# Path to JSBSim files, location of the folders "aircraft", "engines" and "systems"
PATH_TO_JSBSIM_FILES="../../"

fdm = jsbsim.FGFDMExec(PATH_TO_JSBSIM_FILES)  

# Load the F16 aircraft model
fdm.load_model('f16') 

# Set engine running
fdm['propulsion/engine[0]/set-running'] = 1

# State results
results = []

# Initial conditions
fdm['ic/h-sl-ft'] = 20000
fdm['ic/vc-kts'] = 335
fdm['ic/gamma-deg'] = 0

# Initialize the aircraft with initial conditions
fdm.run_ic() 

# Trim
try:
    fdm['simulation/do_simple_trim'] = 1
except RuntimeError as e:
    # The trim cannot succeed. Just make sure that the raised exception
    # is due to the trim failure otherwise rethrow.
    if e.args[0] != 'Trim Failed':
        raise

# Run for 5 min
while fdm.get_sim_time() < 300:
    fdm.run()
    gamma = fdm['flight-path/gamma-deg']
    alt = fdm['position/h-sl-ft']
    pitch = fdm['attitude/theta-deg']
    roll = fdm['attitude/roll-rad'] * 57.29578

    time = fdm.get_sim_time()

    results.append((time, alt, gamma, pitch, roll))

# Extract the results
time, alt, gamma, pitch, roll = zip(*results)

plt.figure()

plt.plot(time, gamma, label='$\gamma$')
plt.plot(time, roll, label='$\phi$')

plt.title("F-16 Trimmed at 20,000ft 335KCAS")
plt.xlabel('Time (s)')
plt.ylabel('Angle (deg)')

plt.legend()
plt.show()

What I noticed while running this is that the root cause doesn't seem to be related to the pitch channel of the FCS, but rather the roll channel.

First off the trim results.

  Full Trim

  Trim successful
  Trim Results:
       Angle of Attack:   1.40  wdot:  1.32e-05 Tolerance: 1e-03  Passed
              Throttle:   0.41  udot: -6.27e-05 Tolerance: 1e-03  Passed
            Pitch Trim:  -0.06  qdot: -6.44e-08 Tolerance: 1e-04  Passed
            Roll Angle:  -0.00  vdot:  2.96e-08 Tolerance: 1e-03  Passed
              Ailerons:  -0.00  pdot:  8.51e-09 Tolerance: 1e-04  Passed
                Rudder:   0.00  rdot: -4.88e-09 Tolerance: 1e-04  Passed

And now the aircraft state (flight path angle, bank angle) over the next 5 minutes.

So even though the trim results set fcs/aileron-cmd-norm to 0, i.e. demanding a zero roll rate command the aircraft starts rolling.

I double-checked the value of fcs/aileron-cmd-norm in case the trim results which only display to 2 decimal places was hiding a very small roll rate command, but at least to 6 decimal places it's 0.

Trim (Event 1) executed at time: 0.108329
    propulsion/engine[0]/n2 = 91.305264
    propulsion/engine[0]/thrust-lbs = 5849.448916
    velocities/vc-kts = 335.000000
    position/h-agl-ft = 20000.000000
    attitude/theta-deg = 1.396447
    attitude/roll-rad = -0.000000
    aero/alpha-deg = 1.396447
    fcs/elevator-cmd-limiter = -0.059655
    fcs/elevator-scheduler = -0.057067
    fcs/pitch-trim-error = -0.017171
    fcs/roll-trim-error = -0.000000
    fcs/aileron-cmd-norm = -0.000000

But by 160 seconds the bank angle has increased to 8.9 deg.

Repeating Notify (Event 3) executed at time: 160.001933
    propulsion/engine[0]/n2 = 91.305264
    propulsion/engine[0]/thrust-lbs = 5853.839479
    velocities/vc-kts = 334.947645
    position/h-agl-ft = 19980.479895
    attitude/theta-deg = 1.280890
    attitude/roll-rad = 0.155827
    aero/alpha-deg = 1.292145
    fcs/elevator-cmd-limiter = -0.059655
    fcs/elevator-scheduler = -0.057260
    fcs/pitch-trim-error = -0.015074
    fcs/roll-trim-error = -0.003813
    fcs/aileron-cmd-norm = 0.000000

The roll channel is a lot simpler than the pitch channel, so it should be easier to debug 😉

[seanmcleod]

As a quick test I changed the proportional gain by an order of magnitude from 3 to 30.

   <pid name="fcs/roll-rate-pid">
     <trigger>fcs/aileron-pid-trigger</trigger>
     <input>fcs/roll-trim-error</input>
     <kp> 30.00000 </kp>
     <ki> 0.00050 </ki>
     <kd> -0.00125 </kd>
   </pid>

With the following results over the 5 minute run.

Repeating Notify (Event 3) executed at time: 300.004666
    propulsion/engine[0]/n2 = 91.305264
    propulsion/engine[0]/thrust-lbs = 5869.799350
    velocities/vc-kts = 336.757080
    position/h-agl-ft = 19886.924149
    attitude/theta-deg = 1.260810
    attitude/roll-rad = 0.000376
    aero/alpha-deg = 1.311653
    fcs/elevator-cmd-limiter = -0.059655
    fcs/elevator-scheduler = -0.057224
    fcs/pitch-trim-error = -0.017606
    fcs/roll-trim-error = -0.000002
    fcs/aileron-cmd-norm = 0.000000

I guess the question is what should the 'optimal' gain be. I'm assuming that users hand flying the F-16 model haven't noticed this since they're fairly active on the stick, and if they do try and fly at a constant bank angle for long periods they probably tweak/correct the bank angle every couple of minutes once it drifts off slightly.

[seanmcleod]

I guess this is a known issue with rate controllers, i.e. no matter how fast your FCS rate controller is, any disturbance e.g. from turbulence will change the state, e.g. initiate a slight roll rate and the FCS will zero out the rate after some non-zero time given the time it takes to process the sensor data and the time it takes to move the actuators etc. But it won't undo any accumulated change in roll angle.

With turbulence I guess on average there will be disturbances in both directions such that they may over time roughly average out to leave you at the same roll angle.

But in the case of JSBSim with no turbulence, which I didn't enable, I'm guessing there is some very small roll rate, possibly even just due to the finite floating point precision and possibly constant in it's direction which is then leading to the results we're seeing.

And I guess typically there is an outer loop to the rate controller, e.g. a pilot or navigation controller etc. which would detect this sort of drift and cancel it out.

[seanmcleod]

The other thing is the PID gains aren't gain scheduled to take into account the different dynamics at different airspeeds/dynamic pressure, so for example this is the result after increasing the airspeed from 335KCAS to 500KCAS. So a much smaller roll excursion, hmm also in the opposite direction...

I think there is some attempt to possibly normalize the roll rate instead of having gain scheduling? Although instead of a fixed gain in this normalization I'd expect to see a reference to dynamic pressure or IAS? @agodemar I'm sure at some point you made mention of this roll rate normalisation being incorrect?

   <!-- Calculate the normalized roll-rate -->
   <pure_gain name="fcs/roll-rate-norm">
    <input>velocities/p-aero-rad_sec</input>
    <gain>0.31821</gain>
   </pure_gain>

[agodemar]

@seanmcleod usually the normalized roll rate is defined as $\hat{p} = pb/(2V)$, and the associated roll coefficient derivative (roll damping derivative) is defined as

$$C_{\mathcal{L}p} = \frac{\partial C{\mathcal{L}}}{\partial \hat{p}} = \frac{\partial C_{\mathcal{L}}}{\partial \left( pb/(2V) \right)} $$

and the contribution to roll moment due to roll rate would be written as
$$\Delta \mathcal{L}(p) = \frac{1}{2}\rho V^2 \ S \ b \ C_{\mathcal{L}_p} \frac{pb}{2V} $$

This means you should have a dependency on airspeed $V$, that is, on Mach number. But maybe the original model includes this effect differently. The above formula is the standard way aerodynamicists report the values of nondimensional derivatives. If the original report on which the model is based complies to the widely accepted definitions, we should check the JSBSim model.

[seanmcleod]

@agodemar taking a look at aero/coefficient/Clp it looks like converting the xml function to a math formula we have this for roll damping?

$$\Delta \mathcal{L}(p) = C_{\mathcal{L}_p}(\alpha) \frac{1}{2}{\rho}V^2sb\frac{b}{2V} p $$

So simplifying to:

$$\Delta \mathcal{L}(p) = C_{\mathcal{L}_p}(\alpha) \frac{1}{2}{\rho}Vs\frac{b^2}{2}p $$

   <function name="aero/coefficient/Clp">
    <description>Roll_moment_due_to_roll_rate_(roll_damping)</description>
    <product>
     <property>aero/qbar-psf</property>
     <property>metrics/Sw-sqft</property>
     <property>metrics/bw-ft</property>
     <property>aero/bi2vel</property>
     <property>velocities/p-aero-rad_sec</property>
       <table>
        <independentVar>aero/alpha-rad</independentVar>
        <tableData>
         -0.1750  -0.3600
         -0.0870  -0.3590
         0.0000  -0.4430
         0.0870  -0.4200
         0.1750  -0.3830
         0.2620  -0.3750
         0.3490  -0.3290
         0.4360  -0.2940
         0.5240  -0.2300
         0.6110  -0.2100
         0.6980  -0.1200
         0.7850  -0.1000
        </tableData>
       </table>
    </product>
   </function>

[seanmcleod]

So taking a look at the roll rate over time, which the roll rate command loop should be keeping at 0 we can see that the roll rate keeps increasing. So there must be a net roll moment even with the roll rate command loop deflecting the ailerons etc. to counter it.

So I looked up all the aerodynamic roll moments that are defined in the F-16's FDM and recorded each moment and calculated the net aerodynamic roll moment.

roll_moments = {
    "aero/coefficient/Clb": [],
    "aero/coefficient/Clb_M": [],
    "aero/coefficient/Clp": [],
    "aero/coefficient/Clr": [],
    "aero/coefficient/Clda": [],
    "aero/coefficient/Clda_M": [],
    "aero/coefficient/Cldr": [],
    "aero/coefficient/Cldr_M": []
}

net_aerodynamic_roll_moment = []

    net_moment = 0
    for moment_property in roll_moments.keys():
        moment = fdm[moment_property]
        net_moment += moment
        roll_moments[moment_property].append(moment)

    net_aerodynamic_roll_moment.append(net_moment)

With the following results.

So the largest moment causing the roll is due to Clb_M, i.e. roll moment due to beta with Mach effect.

   <function name="aero/coefficient/Clb_M">
    <description>Change_in_Roll_moment_due_to_beta_due_to_mach</description>
    <product>
     <property>aero/qbar-psf</property>
     <property>metrics/Sw-sqft</property>
     <property>metrics/bw-ft</property>
     <property>aero/beta-rad</property>
       <table>
        <independentVar>velocities/mach</independentVar>
        <tableData>
         0.6000  0.0000
         0.8000  0.1891
         1.0000  0.1776
         1.2000  0.1375
         1.4000  0.1203
         1.6000  0.0859
        </tableData>
       </table>
    </product>
   </function>

Now what's interesting is if you compare it to Clb in the graph it's of the opposite sign. Is this a typo in the FDM, or is it accurate/correct in terms the roll moment due to beta changing direction/sign past some Mach number? Will need to double-check some of the source documents listed.

Here is Clb from the FDM.

   <function name="aero/coefficient/Clb">
    <description>Roll_moment_due_to_beta</description>
    <product>
     <property>aero/qbar-psf</property>
     <property>metrics/Sw-sqft</property>
     <property>metrics/bw-ft</property>
       <table>
        <independentVar lookup="row">aero/alpha-rad</independentVar>
        <independentVar lookup="column">aero/beta-rad</independentVar>
        <tableData>
                 -0.5240 -0.4360 -0.3490 -0.2620 -0.1750 -0.0870  0.0000  0.0870  0.1750  0.2620  0.3490  0.4360  0.5240
         -0.1750 -0.0090 -0.0070 -0.0000  0.0010  0.0030  0.0010  0.0000 -0.0010 -0.0030 -0.0010  0.0000  0.0070  0.0090
         -0.0870  0.0110  0.0100  0.0100  0.0100  0.0090  0.0040  0.0000 -0.0040 -0.0090 -0.0100 -0.0100 -0.0100 -0.0110
          0.0000  0.0230  0.0230  0.0220  0.0200  0.0170  0.0080  0.0000 -0.0080 -0.0170 -0.0200 -0.0220 -0.0230 -0.0230
          0.0870  0.0370  0.0340  0.0340  0.0300  0.0240  0.0120  0.0000 -0.0120 -0.0240 -0.0300 -0.0340 -0.0340 -0.0370
          0.1750  0.0500  0.0490  0.0470  0.0390  0.0300  0.0160  0.0000 -0.0160 -0.0300 -0.0390 -0.0470 -0.0490 -0.0500
          0.2620  0.0680  0.0630  0.0600  0.0540  0.0410  0.0220  0.0000 -0.0220 -0.0410 -0.0540 -0.0600 -0.0630 -0.0680
          0.3490  0.0890  0.0810  0.0690  0.0570  0.0450  0.0220  0.0000 -0.0220 -0.0450 -0.0570 -0.0690 -0.0810 -0.0890
          0.4360  0.0880  0.0790  0.0670  0.0540  0.0400  0.0210  0.0000 -0.0210 -0.0400 -0.0540 -0.0670 -0.0790 -0.0880
          0.5240  0.0910  0.0600  0.0330  0.0230  0.0160  0.0150  0.0000 -0.0150 -0.0160 -0.0230 -0.0330 -0.0600 -0.0910
          0.6110  0.0760  0.0580  0.0360  0.0060  0.0020  0.0080  0.0000 -0.0080 -0.0020 -0.0060 -0.0360 -0.0580 -0.0760
          0.6980  0.0770  0.0620  0.0350  0.0140  0.0100  0.0130  0.0000 -0.0130 -0.0100 -0.0140 -0.0350 -0.0620 -0.0770
          0.7850  0.0760  0.0590  0.0350  0.0270  0.0190  0.0150  0.0000 -0.0150 -0.0190 -0.0270 -0.0350 -0.0590 -0.0760
        </tableData>
       </table>
    </product>
   </function>

[dpculp]

It looks like the roll command is coming solely from the FCS. Here is a chart showing aileron angle expanded x10 and roll angle. With no other inputs, pilot or turbulence, this is expected.

I tried following the chain of components in the roll FCS, but the mixing of angles and normalized values is throwing me off. It seems the aileron-cmd-norm is used as a proxy for actual roll rate, so we're mixing radians with norms and mixing actual with expected.

[seanmcleod]

My initial suspicion was that there is a small amount of beta developing, which the yaw rate controller doesn't completely dampen, and this induces a roll moment via Clb (roll moment due to beta) and the roll rate controller isn't aggressive enough in damping this roll moment which continues to grow.

But I haven't confirmed it, it could be that the beta is being induced by a small roll command generating a small roll rate, i.e. via Cnp (yaw moment due to roll rate). They both interact/feed each other, but need to figure out which one is the chicken and which one is the egg.

[seanmcleod]

Still not definitive I guess, it could still be the case that a very small roll rate is inducing a beta deviation.

Summary of the rolling moments with Clb and Clb_M combined. Also note the fairly large rolling moment from the rudder, Cldr.

[seanmcleod]

So I decided to comment out the following rolling moments to see what difference it would make.

Clb, Cldr, i.e. rolling moments due to beta and rudder deflection and this is the result over 5min, i.e. not rolling inverted and crashing into the ground.

[dpculp]

Well, that's a vast improvement. Now the focus moves to the yaw channel? Maybe excessive coupling between roll and yaw? I would expect Clb to be pretty high since it's a fighter.

EDIT: The pitch channel is also coupled with the ailerons

[agodemar]

@seanmcleod did you change something in the input files or is this a dissection ot the outputs?

[seanmcleod]

@agodemar for this last graph I took the F-16 FDM and simply commented out the following coefficient 'functions' - Clb and Cldr, i.e. the rolling moments that would be generated due to beta directly i.e. Clb and indirectly Cldr coming from the yaw channel deflecting the rudder to counter beta etc.

[thomasbbrunner]

I've been experimenting with the F-16 model and have also observed similar behavior to what is described here. The airplane seems very prone to rolling, requiring constant input corrections to keep it flying straight and level. This is true also when disabling the fly-by-wire features (i.e., directly actuating the control surfaces).

It also seems that excessive rolling -- such as when performing an aileron roll -- has little to no impact on the aerodynamic drag of the aircraft, which is somewhat counter intuitive.

I'm not very familiar with the aerodynamical concepts discussed here. But are there any takeaways from these experiments? Can we modify the FDM in a way to reduce these issues?

[seanmcleod]

@thomasbbrunner yep, we never really resolved this at the time. I see you're interested in RL, in which case you may be interested to know that for the Darpa Virtual Air Combat Competition - https://github.com/JSBSim-Team/jsbsim/wiki/Interesting-JSBSim-Examples#darpa-virtual-air-combat-competition that JSBSim with the F-16 was used both by the RL teams and also the airforce pilot in the final, without any apparent roll issues 😉

[seanmcleod]

Now the focus moves to the yaw channel?

Yep, so I took a brief look at the yaw FCS channel. First quick test was simply to disconnect the rudder, by renaming the output property.

   <aerosurface_scale name="fcs/rudder-control">
    <input>fcs/rudder-pos-norm</input>
    <range>
     <min>-0.524</min>
     <max>0.524</max>
    </range>
    <output>fcs/rudder-pos-rad-XXX</output>
   </aerosurface_scale>

So with the rudder centered the whole time, I get the same result I did above when I commented out $C_{L_\beta}$ and $C_{L_{\delta_r}}$.

The FCS yaw channel is outputting a small and increasing rudder angle which then generates $\beta$ plus a small rolling moment of it's own. Which given the current roll rate PID controller isn't arrested immediately.

The following shows what is used to calculate the error signal for the FCS yaw channel PID controller.

jsbsim/aircraft/f16/f16.xml

Lines 674 to 702 in 941da03

	<channel name="Yaw">
	<!-- Calculate the normalized yaw-rate -->
	<scheduled_gain name="fcs/yaw-rate-norm">
	<input>velocities/r-aero-rad_sec</input>
	<table>
	<independentVar>velocities/vg-fps</independentVar>
	<tableData>
	80.0 0.0
	100.0 15.0
	150.0 100.0
	</tableData>
	</table>
	</scheduled_gain>
	<!-- Calculate the normalized yaw-load -->
	<pure_gain name="fcs/yaw-load-norm">
	<input>accelerations/n-pilot-y-norm</input>
	<gain>0.25</gain>
	</pure_gain>
	<!--
	- Calculate the difference between the current yaw-rate
	- and the one requiested for.
	-->
	<summer name="fcs/yaw-trim-error">
	<input>fcs/rudder-cmd-norm</input>
	<input>fcs/yaw-rate-norm</input>
	<input>fcs/yaw-load-norm</input>
	</summer>

Having a quick glance at some F-16 documentation elsewhere I'm not sure this is correct at a high-level, i.e. from what I've read so far the pilot doesn't command a yaw rate from the rudder pedals.

https://arc.aiaa.org/doi/book/10.2514/4.867873

https://ntrs.nasa.gov/api/citations/19800005879/downloads/19800005879.pdf

https://www.falcon-bms.com/wp-content/uploads/2021/08/FM_Developers_Notes_Part_4.pdf