Code covered by the BSD License

# H-infinity sub-optimal causal FIR inverse filter via LMI

### Masaaki Nagahara (view profile)

This computes the optimal inverse filter minimizing the H-infinity norm of the error system via LMI.

dfdesign_w_lmi(phi, w, d, n);
```function [psi gopt] = dfdesign_w_lmi(phi, w, d, n);

% [psi gopt] = dfdesign(phi, d, n);
% DFDESIGN computes the H-infinity optimal inverse FIR filter of phi(z).
% The resulting filter minimizes the H-infinity norm of the error system
%     E_w(z) = [z^(-d) - psi(z)phi(z)]w(z).
% The LMI method based on KYP lemma is used in this function.
%
% [INPUT]
% phi: target system (discrete-time)
% w: weighting function
% d: delay
% n: length of the FIR filter
%
% [OUTPUT]
% psi: resulting inverse FIR filter (in state-space representation)
% gopt: optimal value

Ts = phi.Ts;
delayden = zeros(1, d+1);
delayden(1) = 1;
delay = tf(1, delayden, Ts);

[Aphi Bphi Cphi Dphi] = ssdata(ss(phi, 'min')*w);
[Ad Bd Cd Dd] = ssdata(ss(delay, 'min')*w);

Apsi11 = zeros(n-1,1);
Apsi21 = 0;
Apsi12 = eye(n-1);
Apsi22 = zeros(1,n-1);
Apsi = [Apsi11 Apsi12; Apsi21 Apsi22];
Bpsi = zeros(n, 1);
Bpsi(n) = 1;
Kpsi1 = eye(n);
Kpsi2 = zeros(1,n);
Kpsi = [Kpsi1; Kpsi2];
Mpsi = zeros(n+1,1);
Mpsi(size(Mpsi,1)) = 1;

A13 = zeros(size(Apsi, 1), size(Ad, 2));
A21 = zeros(size(Aphi, 1), size(Apsi, 2));
A23 = zeros(size(Aphi, 1), size(Ad, 2));
A31 = zeros(size(Ad, 1), size(Apsi, 2));
A32 = zeros(size(Ad, 1), size(Aphi, 2));
A = [Apsi Bpsi*Cphi A13; A21 Aphi A23; A31 A32 Ad];
B = [Bpsi*Dphi; Bphi; -Bd];
Cp13 = zeros(size(Kpsi, 1), size(Cd, 2));
Cp = [Kpsi Mpsi*Cphi Cp13];
Cq1 = zeros(1, size(Kpsi, 2));
Cq2 = zeros(1, size(Mpsi*Cphi, 2));
Cq = [Cq1 Cq2 Cd];
Dp = Mpsi*Dphi;
I22 = eye(size(B, 2), size(B, 2));
I33 = eye(1, 1);
Pdim = size(A, 2);

setlmis([]);
a=lmivar(2,[1 n+1]);
P=lmivar(1,[Pdim 1]);
x=lmivar(2,[1 1]);

lmiterm([1 1 1 P],A',A);                        % LMI #1: A'*P*A
lmiterm([1 1 1 P],1,-1);                        % LMI #1: -P
lmiterm([1 2 1 P],B',A);                        % LMI #1: B'*P*A
lmiterm([1 2 2 x],.5*1,-I22,'s');               % LMI #1: -x*I22 (NON SYMMETRIC?)
lmiterm([1 2 2 P],B',B);                        % LMI #1: B'*P*B
lmiterm([1 3 1 a],1,Cp);                        % LMI #1: a*Cp
lmiterm([1 3 1 0],Cq);                          % LMI #1: Cq
lmiterm([1 3 2 a],1,Dp);                        % LMI #1: a*Dp
lmiterm([1 3 3 x],.5*1,-I33,'s');               % LMI #1: -x*I33 (NON SYMMETRIC?)

lmiterm([-2 1 1 0],P);                          % LMI #2: P

lmisys=getlmis;

% mincx
lmin = decnbr(lmisys);
lmic = zeros(lmin, 1);
for j=1:lmin,
[xj] = defcx(lmisys,j,x);
lmic(j) = xj;
end
[gopt, xopt] = mincx(lmisys, lmic);

aopt = dec2mat(lmisys, xopt, 1);
Cpsi = aopt * Kpsi;
Dpsi = aopt * Mpsi;
psi = ss(Apsi, Bpsi, Cpsi, Dpsi, Ts);

```