system identification: computing output of a model

Dear All,

I use system identification toolbox to estimate parameters of an ARX model and get, for example:

A(q) = 1 - 1.5 q^-1 + 0.7 q^-2

B(q) = q^-1 + 0.5 q^-2 - 0.1 q^-3

which equal to the following difference equation:

y(k) - 1.5 y(k-1) + 0.7 y(k-2) = u(k-1) + 0.5 u(k-2) - 0.1 u(k-3)

To compute output at time k, let say k = 4, I need to supply y(k-1),y(k-2)u(k-1),u(k-2), and u(k-3).

The question is if I use function 'sim' to compute output of the system above, how can set the initial value? Until now, I cannot make the Matlab's output is the same as hand calculation.

Furthermore, if I converted the model into state space form, it seems that the input matrix B from Matlab is not the same as those of standard signals and systems textbooks when I converted it into the canonic observable form.

Probably, I am missing something here.

Thank you for your help in advance.

Zul

Answer 1

[Rajiv Singh]

Specifying initial conditions for simulation of polynomial models

The SIM command does not support this directly. But there are two options to do this:

1. Simulate using FILTER command
2. Simulate by converting to state-space form first

Using filter: See the documentation of FILTER; Use B = model.B; A = model.F and generate initial conditions based on the assumed Direct Form II Transposed filter structure.

Using state space:

1. Convert the model into state-space form using IDSS or SS command.
2. Map the input-output data history to the initial states x0 of the state-space model obtained above. Use the function DATA2STATE for this (see below).
3. Use SIM command on the state-space model with the initial conditions x0.

Form of state space model You should not rely on the form of the state-space model returned on conversion from a transfer function or polynomial model. There is no unique representation. Some popular forms are facilitated by the CANON command in Control System Toolbox that might help you convert your model into a textbook canonical form. Use the "SS" command to convert an IDPOLY model into the state-space model of Control System Toolbox: sys = ss(idpoly_model(:,'measured'))

% CODE START
function x0 = data2state(sys, Init)
% Map data history to initial states of discrete-time state-space model. 
% Syntax:
%     X0 = DATA2STATE(SYS, INIT)
%
% SYS: Discrete-time state-space model with no internal delays. 
% Init: Structure with fields 'Input' and 'Output'. If Input and/or Output
%       contains single row of values, scalar expansion is used; the supplied
%       signal values are assumed to be constant over time. For exact
%       match, Init.Input and Init.Output must contain as many samples as
%       there are states in sys. Otherwise, least-squares (if #samples>nx)
%       or minimum-norm solution (if #samples<nx) is obtained. INIT.Input
%       and INIT.Output data must correspond to time instants (-Nx, -Nx+1,
%       ..., -1), the last sample being the latest.
%
% X0: Value of states at the sample after the time instant of the last data
% sample in INIT. Column vector.
% 
% Copyright 2011 The MathWorks, Inc.

% data preparation
assert(sys.Ts~=0,'Need discrete-time model');
[A,B,C,D] = ssdata(sys);
[ny, nx] = size(C);  nu = size(B,2); 
y = Init.Output; u = Init.Input;
if isempty(y), y = zeros(1,ny); end
if isempty(u), u = zeros(1,nu); end
if isrow(y) && nx>1, y = repmat(y,[nx,1]); end
if isrow(u) && nx>1, u = repmat(u,[nx,1]); end

% initial state computation
x00 = x0est([y,u],A,B,C,D,zeros(nx,ny),ny,nu,250e3,eye(ny));
u = [u; zeros(1,nu)];
x = ltitr(A,B,u,x00);
x0 = x(end,:).';

% CODE END