Re: control problem: lsim

[Lukas Reichlin]

On 25.02.2013, at 13:36, Rudolf Widmer-Schnidrig <address@hidden> wrote:

> On 25.02.13 07:45, Lukas Reichlin wrote:
>> On 24.02.2013, at 23:17, Rudolf Widmer-Schnidrig <address@hidden> wrote:
>> 
>>> Dear Lukas, dear List,
>>> 
>>> What is the best way to evaluate the quality of a LTI model
>>> which I have obtained like this
>>> 
>>> octave> dat = iddata(yout,xin);      % yout is observed output to input xin
>>> octave> sys = arx (dat,'NA',2,'NB',1')
>>> 
>>> From the control manual I expected that I could run
>>> 
>>> octave> lsim(sys,xin)
>>> 
>>> to predict the output for input xin given the system sys.
>>> 
>>> but I get the error: lsim: input vector u must have 2 columns
>>> 
>>> what am  missing?
>>> 
>>>       many thanks in advance,             -Ruedi
>> Hi Ruedi,
>> 
>> If dat has p outputs and m inputs,
>> 
>> [n, p, m, e] = size (dat)
>> [n, p] = size (yout)
>> [n, m] = size (xin)
>> 
>> sys becomes an p by m+p LTI model (p outputs and m+p inputs).
>> The first m inputs (1:m) are u(t) and the remaining p inputs (m+(1:p)) are 
>> e(t)
>> 
>>           A(q) y(t) = B(q) u(t) + e(t)
>> 
>> You are only interested in u(t), therefore select the inputs by  sys(:, 1:m)
>> 
>> lsim (sys(:, 1:m), xin)
>> 
>> Maybe I should improve the documentation of arx.
>> 
>> Best regards,
>> Lukas
>> 
>> 
>> 
> Dear Lukas,
> 
> thank you for the clarifying words!
> 
> I deal with the simplest case: single input and single output. Thus xin and 
> yout are plain vectors.
> 
> With:
> octave> [yA,tA,xA] = lsim(sys(:,1),xin);
> 
> I get the simulated output of my system: yA    (is this interpretation of the 
> output correct?)
> 
> and by plotting
> 
> octave> plot(ta,yout,'k;measured output;',ta,yA,'r;simulated output;')
> 
> I can see how well the LTI system (sys) is able to mimic the physical world. 
> (actually not too bad for the beginning)
> 
> Now I want to convert the discrete time model to a continuous time model with 
> "d2c".
> 
> I start with specifying the sample interval dt=0.01s
> 
> octave> set(sys,"tsam",0.01);
> 
> then do the conversion
> 
> octave> Csys = d2c(sys,'bilin');
> 
> and run into the next problem if I want to see the contents of Csys:
> 
> octave> octave:44> Csys
> 
> error: lti: filtdata: require discrete-time system
> error: called from:
> error:   /Users/widmer/octave/control-2.4.0/@lti/filtdata.m at line 62, 
> column 5
> error:   /Users/widmer/octave/control-2.4.0/@tf/display.m at line 36, column 
> 17
> 
> while: octave> get(Csys,"num") returns sensible values.

There was a bug in the tf display routine. It is fixed in control-2.4.2. Please 
update your package.

> 
> 
> Then more general questions:
> The best model  I have to describe my system is H(s)=G*s/((s-p1)(s-p2))   
> where p1 and p2 are a pair of complex conjugate poles and G is a real number.
> So far the models I get from arx all look like low-pass filters while my 
> system is a band pass (as is obvious from a spectral division of 
> output/input).
> To evaluate the frequency response of sys I did:
> 
> octave> sys = arx (dat,'NA',2,'NB',1);
> set(sys,"tsam",0.01);
> [H,w,tsam]=frdata(sys(:,1),"v");
> loglog(w,abs(H))
> 
> Am I not setting NA and NB correctly?
> 
> Is there a way to force arx to put the zeros at s=0
> 
> Is there a way to estimate the state-space model which is closest to some 
> starting model?
> 
>           many thanks             -Ruedi

Does this help?

        [sys, x0] = arx (dat, …)

produces a state-space model sys and initial state vector x0. Also see Thomas' 
mail.

Lukas