Linear Approximation of Complex Systems by Identification

This example shows how to obtain linear approximations of a complex, nonlinear system by means of linear model identification. The approach is based on selection of an input signal that excites the system. A linear approximation is obtained by fitting a linear model to the simulated response of the nonlinear model for the chosen input signal.

Introduction

In many situations, a linear model is obtained by simplification of a more complex nonlinear system under certain local conditions. For example, a high fidelity model of the aircraft dynamics may be described by a detailed Simulink model. A common approach taken to speed up the simulation of such systems, study their local behavior about an operating point or design compensators is to linearize them. If we perform the linearization of the original model analytically about an operating point, this, in general, would yield a model of an order as high as, or close to, the number of states in the original model. This order may be unnecessarily high for the class of inputs it is supposed to be used with for analysis or control system design. Hence we can consider an alternative approach centered around collecting input-output data from the simulation of the system and using it to derive a linear model of just the right order.

Analytical Linearization of F14 Model

Consider the F14 model. This is a already a linear model but contains derivative blocks and sources of disturbance that may affect the nature of its output. We can "linearize" it between its one input port and the two output ports as follows:

```open_system('idF14Model') IO = getlinio('idF14Model') syslin = linearize('idF14Model',IO) ```
```3x1 vector of Linearization IOs: -------------------------- 1. Linearization input perturbation located at the following signal: - Block: idF14Model/Pilot - Port: 1 - Signal Name: Stick Input 2. Linearization output measurement located at the following signal: - Block: idF14Model/Gain5 - Port: 1 - Signal Name: Angle of Attack 3. Linearization output measurement located at the following signal: - Block: idF14Model/Pilot G force - Port: 1 - Signal Name: Pilot G force syslin = A = Transfer Fcn Derivative Transfer Fcn Derivative1 Transfer Fcn -0.6385 0 689.4 0 Derivative 1 -1e+05 0 0 Transfer Fcn -0.00592 0 -0.6571 0 Derivative1 0 0 1 -1e+05 Actuator Mod 0 0 1.424 0 Alpha-sensor 0.001451 0 0 0 Stick Prefil 0 0 0 0 Pitch Rate L 0 0 1 0 Proportional 0 0 -0.8156 0 Actuator Mod Alpha-sensor Stick Prefil Pitch Rate L Transfer Fcn -1280 0 0 0 Derivative 0 0 0 0 Transfer Fcn -137.7 0 0 0 Derivative1 0 0 0 0 Actuator Mod -20 2.986 -39.32 -1.67 Alpha-sensor 0 -2.526 0 0 Stick Prefil 0 0 -22.52 0 Pitch Rate L 0 0 0 -4.144 Proportional 0 -1.71 22.52 0.9567 Proportional Transfer Fcn 0 Derivative 0 Transfer Fcn 0 Derivative1 0 Actuator Mod -3.864 Alpha-sensor 0 Stick Prefil 0 Pitch Rate L 0 Proportional 0 B = Stick Input Transfer Fcn 0 Derivative 0 Transfer Fcn 0 Derivative1 0 Actuator Mod 0 Alpha-sensor 0 Stick Prefil 1 Pitch Rate L 0 Proportional 0 C = Transfer Fcn Derivative Transfer Fcn Derivative1 Angle of Att 0.001451 0 0 0 Pilot G forc -3106 3.106e+08 7.083e+04 -7.081e+09 Actuator Mod Alpha-sensor Stick Prefil Pitch Rate L Angle of Att 0 0 0 0 Pilot G forc 0 0 0 0 Proportional Angle of Att 0 Pilot G forc 0 D = Stick Input Angle of Att 0 Pilot G forc 0 Continuous-time state-space model. ```

syslin is a model with 2 outputs, one input and 9 states. This is because the linearization path from the "Stick Input" input to the two outputs has 9 states in the original system. We can verify this by using `operpoint`:

```operpoint('idF14Model') ```
```ans = Operating point for the Model idF14Model. (Time-Varying Components Evaluated at time t=0) States: ---------- x _ (1.) idF14Model/Actuator Model 0 (2.) idF14Model/Aircraft Dynamics Model/Transfer Fcn.1 0 (3.) idF14Model/Aircraft Dynamics Model/Transfer Fcn.2 0 (4.) idF14Model/Controller/Alpha-sensor Low-pass Filter 0 (5.) idF14Model/Controller/Pitch Rate Lead Filter 0 (6.) idF14Model/Controller/Proportional plus integral compensator 0 (7.) idF14Model/Controller/Stick Prefilter 0 (8.) idF14Model/Dryden Wind Gust Models/Q-gust model 0 (9.) idF14Model/Dryden Wind Gust Models/W-gust model 0 0 Inputs: None ---------- ```

Could this order be reduced while still maintaining fidelity to the responses for the chosen square-wave ("Stick input") input?

Preparing Identification Data

Simulate the model and log the input square wave (u) and the outputs "Angle of attack" (y1) and "Pilot G force" (y2) for 0:30 second time span. This data, after interpolation to a uniformly spaced time vector (sample time of 0.0444 seconds), is stored in the "idF14SimData.mat" file.

```load idF14SimData Z = iddata([y1, y2],u,Ts,'Tstart',0); Z.InputName = 'Stick input'; Z.InputUnit = 'rad/s'; Z.OutputName = {'Angle of attack', 'Pilot G force'}; Z.OutputUnit = {'rad', 'g'}; t = Z.SamplingInstants; subplot(311) plot(t,Z.y(:,1)), ylabel('Angle of attack (rad)') title('Logged Input-Output Data') subplot(312) plot(t,Z.y(:,2)), ylabel('Pilot G force (g)') subplot(313) plot(t,Z.u), ylabel('Stick input (rad/s)') axis([0 30 -1.2 1.2]) xlabel('Time (seconds)') ```

Estimation of State Space Models

Estimate state-space models of orders 2 to 4 using the `ssest` command. We configure the estimation to use "simulation" focus and choose to not estimate the disturbance component of the model.

```opt = ssestOptions('Focus','simulation'); syslin2 = ssest(Z, 2, 'DisturbanceModel', 'none', opt); syslin3 = ssest(Z, 3, 'DisturbanceModel', 'none', opt); syslin4 = ssest(Z, 4, 'DisturbanceModel', 'none', opt); ```

Compare the performance of the linearized model `syslin` and the three identified models to the data. Note that syslin is an SS model while `syslin2`, `syslin3` and `syslin4` are IDSS models.

```syslin.InputName = Z.InputName; syslin.OutputName = Z.OutputName; % reconcile names to facilitate comparison clf compare(Z, syslin, syslin2, syslin3, syslin4) ```

The plot shows that the 3rd order model (`syslin3`) works pretty well as a linear approximation of aircraft dynamic about its default (t=0) operating conditions. It fits the data slightly better than the one obtained by analytical linearization (`syslin`). If the original idF14Model is linear, why doesn't its linearization result, `syslin`, produce a 100% fit to the data? There are two reasons for this:

1. The measured outputs are affected by the wind gusts which means that the logged outputs are not simply a function of the Stick input. There are disturbances affecting it.

2. The "Pilot G force" block uses derivative blocks whose linearization depends upon the value of the time constant "c". "c" should be small (we use 1e-5) but not zero. A non-zero value for c introduces an approximation error in the linearization.

Let us view the parameters of the model `syslin3` which seems to capture the responses pretty well:

```syslin3 ```
```syslin3 = Continuous-time identified state-space model: dx/dt = A x(t) + B u(t) + K e(t) y(t) = C x(t) + D u(t) + e(t) A = x1 x2 x3 x1 -1.006 2.029 0.5842 x2 -8.284 -19.39 5.611 x3 2.784 12.63 -6.956 B = Stick input x1 0.2614 x2 -5.512 x3 3.606 C = x1 x2 x3 Angle of att -8.841 0.5347 1.402 Pilot G forc -86.42 15.85 66.12 D = Stick input Angle of att 0 Pilot G forc 0 K = Angle of att Pilot G forc x1 0 0 x2 0 0 x3 0 0 Parameterization: FREE form (all coefficients in A, B, C free). Feedthrough: none Disturbance component: none Number of free coefficients: 18 Use "idssdata", "getpvec", "getcov" for parameters and their uncertainties. Status: Estimated using SSEST on time domain data "Z". Fit to estimation data: [98.4;97.02]% FPE: 2.367e-05, MSE: 0.1103 ```

Model Simplification by Reduction and Estimation

We can also take the approach of reducing the order of the linearized model `syslin` and refining the parameters of the reduced model to best fit the data Z. To figure out a good value for the reduced order, we use `hsvd`:

```[S, BalData] = hsvd(syslin); clf; bar(S) ```

The bar chart shows that the singular values are quite small for states 4 and above. Hence 3rd order might be optimal for reduction.

```sysr = balred(syslin,3,BalData) opt2 = bodeoptions; opt2.PhaseMatching = 'on'; clf; bodeplot(sysr,syslin,opt2) ```
```sysr = A = x1 x2 x3 x1 -2.854 -7.61 -54.04 x2 -0.9714 2.341 9.123 x3 0.6979 -7.203 -24.08 B = Stick input x1 -137.7 x2 -869 x3 -506.7 C = x1 x2 x3 Angle of att -0.0005063 -0.0008826 -0.001016 Pilot G forc -0.005926 -0.04692 -0.1646 D = Stick input Angle of att -0.03784 Pilot G forc -1.617 Continuous-time state-space model. ```

The bode plot shows good fidelity up to 10 rad/s. `sysr` is able to emulate the responses as good as the original 9-state model as shown by the `compare` plot:

```compare(Z, sysr, syslin) ```

Let us refine the parameters of `sysr` to improve its fit to data. For this estimation, we choose the "Levenberg-Marquardt" search method and change the maximum number of allowable iterations to 10. These choices were made based on some trial and error. We also turn on the display of estimation progress.

```opt.Display = 'on'; opt.SearchMethod = 'lm'; opt.SearchOptions.MaxIterations = 10; sysr2 = ssest(Z, sysr, opt) compare(Z, sysr2) ```
```sysr2 = Continuous-time identified state-space model: dx/dt = A x(t) + B u(t) + K e(t) y(t) = C x(t) + D u(t) + e(t) A = x1 x2 x3 x1 -4.048 -7.681 -54.01 x2 -0.4844 1.549 8.895 x3 -0.2398 -6.777 -25.78 B = Stick input x1 -137.7 x2 -869 x3 -506.7 C = x1 x2 x3 Angle of att -0.0003361 -0.0004964 0.00215 Pilot G forc -0.01191 -0.03599 -0.1434 D = Stick input Angle of att 0.003022 Pilot G forc 0.6438 K = Angle of att Pilot G forc x1 0 0 x2 0 0 x3 0 0 Parameterization: FREE form (all coefficients in A, B, C free). Feedthrough: yes Disturbance component: none Number of free coefficients: 20 Use "idssdata", "getpvec", "getcov" for parameters and their uncertainties. Status: Estimated using SSEST on time domain data "Z". Fit to estimation data: [98.78;97.03]% FPE: 1.434e-05, MSE: 0.1097 ```

The refined model `sysr2` fits the response of the F14 model quite well (about 99% for the first output, and about 97% for the second).

Conclusions

We showed an alternative approach to analytical linearization for obtaining linear approximations of complex systems. The results are derived for a particular input signal and, strictly speaking, are applicable to that input only. To improve the applicability of results to various input profiles, we could perform several simulations using the various types of inputs. We could then merge the resulting datasets into one multi-experiment dataset (see `iddata/merge`) and use it for estimations. In this example, we used a complex, but linear system for convenience. The real benefit of this approach would be seen for nonlinear systems.

We also showed an approach for reducing the order of a linear system while keeping the reduced model faithful to the simulated response of the original Simulink model. The role of the reduced model `sysr` was to provide an initial guess for the estimated model `sysr2`. The approach also highlights the fact that any linear system, including those of a different class, can be used as the initial model for estimations.

```bdclose('idF14Model') ```