## Example 17: Discretization of continuous LPV model of F-16

```close all; clear; clc;
```

## F-16 VISTA Flight Model

For model see [1]

```% continuous-model
A0 = [0.22 1; 17.1 -0.228];
A1 = [-4.10e-7 0; -8.07e-4 7.06e-6];
A2 = [-2.6 0; -68.4 -2.12];
A3 = [5.15e-5 0; 3.31e-3 4.86e-5];
A4 = [0 0; 56.2 0];
A5 = [0 0; -2.92e-3 0];
B0 = [-1.38e-3; -8.16];
B1 = [8.75e-8; 1.73e-4];
B2 = [-0.34; 40.6];
B3 = [7.89e-6; -8.96e-4];
B4 = [0; -99.3];
B5 = [0; 2.42e-3];
C0 = eye(2);

% make affine lpv system
sys = idafflpv([A0 A1 A2 A3 A4 A5],[B0 B1 B2 B3 B4 B5],[C0 zeros(2,5*2)],zeros(2,6),...
'statename',{'alpha' 'q'},'schedulingname',{'h' 'M' '(h*M)' '(M^2)' '(h*M^2)'},...
'inputname','elevator','outputname',{'angle-of-attack' 'pitch-rate'},'Name','F16 Aircraft');
```

## Discretization of continuous LPV model

```% discrete approximation
sysde1 = c2d(sys,0.05,'euler','linear');
sysde3 = c2d(sys,0.05,'tustin','linear');
```

Discretization results

```% scheduling signals
t = (0:0.05:1)';
M = linspace(0.4,0.6,length(t))';
h = linspace(5000,25000,length(t))';
p = [h M h.*M M.^2 M.^2.*h];
pq = [h M h.*M M.^2 M.^2.*h  h.*h h.*M h.*h.*M h.*M.^2 h.*M.^2.*h...
M.*M M.*h.*M M.*M.^2 M.*M.^2.*h  h.*M.*h.*M h.*M.*M.^2 h.*M.*M.^2.*h...
M.^2.*M.^2 M.^2.*M.^2.*h  M.^2.*h.*M.^2.*h];

% pole and zero locations
sys0 = idafflpv2ss(sys,p);
figure, pzmap(sys0);
syse1 = idafflpv2ss(sysde1,p);
syse2 = idafflpv2ss(sysde2,pq);
syse3 = idafflpv2ss(sysde3,p);
syse4 = idafflpv2ss(sysde4,pq);
figure, pzmap(syse1,'y',syse2,'g',syse3,'c',syse4,'b');
```

Performance of discretization

```% sample times
Ts = [0.01 0.02 0.05 0.07 0.1 0.2 0.5 0.7 1];

% scheduling signals
t = 0:Ts(1):100;
u = sin(0.2.*t)';
M = linspace(0.4,0.6,length(t))';
h = linspace(5000,25000,length(t))';
p = [h M h.*M M.^2 M.^2.*h];
pq = [h M h.*M M.^2 M.^2.*h  h.*h h.*M h.*h.*M h.*M.^2 h.*M.^2.*h...
M.*M M.*h.*M M.*M.^2 M.*M.^2.*h  h.*M.*h.*M h.*M.*M.^2 h.*M.*M.^2.*h...
M.^2.*M.^2 M.^2.*M.^2.*h  M.^2.*h.*M.^2.*h];

% continuous simulation
[yc,tc] = sim(sys,u,t,p);

% euler approximation
sysde1 = c2d(sys,Ts(1),'euler','linear');
[ye1,te1] = sim(sysde1,u,t,p);
[yeq1,teq1] = sim(sysdeq1,u,t,pq);

n = 1;
VAFZ = zeros(2,length(Ts));
VAFE = zeros(2,length(Ts));
VAFT = zeros(2,length(Ts));
VAFZ1 = zeros(2,length(Ts));
VAFE1 = zeros(2,length(Ts));
VAFT1 = zeros(2,length(Ts));
VAFEQ = zeros(2,length(Ts));
VAFTQ = zeros(2,length(Ts));
VAFE1Q = zeros(2,length(Ts));
VAFT1Q = zeros(2,length(Ts));
for j = 1:length(Ts);
t = 0:Ts(j):100;
u = sin(0.2.*t)';
M = linspace(0.4,0.6,length(t))';
h = linspace(5000,25000,length(t))';
p = [h M h.*M M.^2 M.^2.*h];
pq = [h M h.*M M.^2 M.^2.*h  h.*h h.*M h.*h.*M h.*M.^2 h.*M.^2.*h...
M.*M M.*h.*M M.*M.^2 M.*M.^2.*h  h.*M.*h.*M h.*M.*M.^2 h.*M.*M.^2.*h...
M.^2.*M.^2 M.^2.*M.^2.*h  M.^2.*h.*M.^2.*h];

% symbolic
sysdz = c2d(sys,Ts(j),'zoh','linear');
[yz,tz] = sim(sysdz,u,t',p);

% euler approximation
sysde = c2d(sys,Ts(j),'euler','linear');
[ye,te] = sim(sysde,u,t',p);

% tustin approximation
sysdt = c2d(sys,Ts(j),'tustin','linear');
[yt,tt] = sim(sysdt,u,t',p);

% euler approximation
[yeq,teq] = sim(sysdeq,u,t',pq);

% tustin approximation
[ytq,ttq] = sim(sysdtq,u,t',pq);

ys = interp1q(tc,yc',t');
ys1 = interp1q(te1,ye1,t');
ysq1 = interp1q(teq1,yeq1,t');
VAFZ(:,j) = vaf(ys,yz);
VAFE(:,j) = vaf(ys,ye);
VAFT(:,j) = vaf(ys,yt);
VAFZ1(:,j) = vaf(ys1,yz);
VAFE1(:,j) = vaf(ys1,ye);
VAFT1(:,j) = vaf(ys1,yt);
VAFEQ(:,j) = vaf(ys,yeq);
VAFTQ(:,j) = vaf(ys,ytq);
VAFE1Q(:,j) = vaf(ysq1,yeq);
VAFT1Q(:,j) = vaf(ysq1,ytq);
end

figure,semilogx(Ts,VAFEQ(1,:),'y',Ts,mean(VAFE),'g',Ts,mean(VAFT),'c',Ts,mean(VAFTQ),'b',Ts,mean(VAFZ),'r');
xlabel('Sample Time [s]')
ylabel('VAF [%]')

figure,semilogx(Ts,mean(VAFE1),'y',Ts,mean(VAFE1Q),'g',Ts,mean(VAFT1),'c',Ts,mean(VAFT1Q),'b',Ts,mean(VAFZ1),'r');
xlabel('Sample Time [s]')
ylabel('VAF [%]')

figure,semilogx(Ts,mean(VAFE),'y',Ts,mean(VAFEQ),'g',Ts,mean(VAFT),'c',Ts,mean(VAFTQ),'b',Ts,mean(VAFZ),'r');
xlabel('Sample Time [s]')
ylabel('VAF [%]')