immense difference between continous and discrete state space model

Hey guys,

I'm trying to set up a dc motor position control by using a full state-feedback controller with integral action. As long as i stick to a time continuous design everything is fine, but as soon as I transform the models into a discrete form the step responses go wild. The attached source code just creates the matrices and gains that are used in a simulink model. The step command gives an impression of what I mean.

Ra  = 2.06;            
La = 0.000238;         
k_phi = 1/((406*2*pi)/60); 
J = 1.07e-6;           
Bm = 12e-7;             

Ts=100e-4;

%general model
%states: theta, dtheta, i, derr
A = [0 1 0;
    0 -Bm/J k_phi/J;
    0 -k_phi/La -Ra/La];
b = [0 ; 0 ; 1/La];
cT = [1 0 0];
d = [0];
asys=ss(A,b,cT,d);           
dsys=c2d(asys,Ts,'zoh');     

%model for integral action
Ai = [0 1 0 0;
    0 -Bm/J k_phi/J 0;
    0 -k_phi/La -Ra/La 0;
    1 0 0 0];
bi = [0 ; 0 ; 1/La ; 0 ];
br = [0 ; 0 ; 0; -1];
cTi = [1 0 0 0];
di = [0];
asysi=ss(Ai,bi,cTi,di);      
dsysi=c2d(asysi,Ts,'zoh');   

% feedback controller
p1 = -100;
p2 = -110;
p3 = -200;
p4 = -300;
poles_k = [p1 p2 p3 p4];                
poles_k = exp(Ts*poles_k);        
K = place(dsysi.a,dsysi.b,poles_k);    
Ki=K(4);                                
K=K(1:3);                               

% observer
poles_obsv = 1000*poles_k(1:3);         
poles_obsv = exp(Ts*poles_obsv);   
L=place(dsys.a',dsys.c',poles_obsv);   
L=L';

t = 0:0.001:.1;
sys_cl = ss(Ai-bi*[K Ki],br,cTi,di);
step(sys_cl,t)
hold on;
sys_cl = ss(dsysi.a-dsysi.b*[K Ki],br,dsysi.c,dsysi.d);
step(sys_cl,t)

All the matrices of dsys are implemented in the simulink model. If one takes the matrices of asys instead the system behaviour is sufficient... Has anybody a clue of what is going wrong here?

br Kev

Answer 1

[Azzi Abdelmalek]

If you can simulate your two systems using Simulink it will be easier to detect the problem. What you did, is find controller parameter in continuous domain, using discrets poles which are unstable, for continuous modeles

poles_k = [p1 p2 p3 p4];                
poles_k = exp(Ts*poles_k);        
K = place(dsysi.a,dsysi.b,poles_k);

poles_k =

    1.1052    1.1052    1.1052    1.1052 %positives unstables poles

% you simulate your continuous model with above unstables parameters,

 sys_cl = ss(Ai-bi*[K Ki],br,cTi,di);

 step(10*sys_cl,t)   % even it works, it's sensless

I have just noticed that your initial model is unstable, which is wierd for a dc machine. I suggest then to post your simulink model

In your code, the observer is not used

[Azzi Abdelmalek]

other remark, in your simulation, there is no any feedback. It's better if you separate your initial system from the controller, it's easy to do it with simulink, you can also do it with matlab by using append and other command

[Babak]

Ai-bi*[K Ki] is the feedback state space matrix of Ai

[Azzi Abdelmalek]

no it's not. It's just an equivalent model, and there is no any feedback

[Babak]

sys_cl = ss(Ai-bi*[K Ki],br,cTi,di)

is the closed loop system with a feedback loop of the open loop system of

sys_cl = ss(Ai,br,cTi,di)

Using K = place(.), the pole placemnt method is used to design the closed loop poles.

[Azzi Abdelmalek]

1. global system(Initial system+Integral action ) is not correct
2. the pôles used to obtain gain K are unstables(>0)
3. the initial system is unstable, which make it difficult to control
4. when you use step(sys_cl,t). there is only one input and one output, how can you imagine there is any feedback. To simulate a feedback, you have to define two blocks: one for initial system and another for the controller u(t)=-k1 x1- k2 x2-k3 x3-ki xi. and create links between them and yref by using append and connect commands.

[Kevin]

the given poles are all placed in the s-plane:

poles_k = [p1 p2 p3 p4];

now they are transformed into the z-plane for discretisation

poles_k = exp(Ts*poles_k);        
K = place(dsysi.a,dsysi.b,poles_k);

considering the z-plane the poles are all within the unit circle, so they should be stable. the step response was just meant to compare the continuous and the discrete approaches for integral action.

Answer 2

[Babak]

Do the responses of the open loop systems look similar?

Just for your information, MATLAB's discrete integral and derivative gains are not equivallent to the corresponding gains of the continuous system. You need to multiply or divide them by Ts. You may want to try this too.

[Azzi Abdelmalek]

c2d is handling the discretization

Answer 3

[Babak]

I edited your code:

 Ra  = 2.06;            
 La = 0.000238;         
 k_phi = 1/((406*2*pi)/60); 
 J = 1.07e-6;           
 Bm = 12e-7;             
 Ts=100e-4;

 %model for integral action
Ai = [0 1 0 0;
    0 -Bm/J k_phi/J 0;
    0 -k_phi/La -Ra/La 0;
    1 0 0 0];
 bi = [0 ; 0 ; 1/La ; 0 ];
 br = [0 ; 0 ; 0; -1];
 cTi = [1 0 0 0];
 di = [0];
 asysi=ss(Ai,bi,cTi,di)    ;  
 dsysi=c2d(asysi,Ts,'zoh') ;  
 % p1 = pole(asysi)
 % exp(Ts.*p1)
 % pole(dsysi)
 % return
 % step(asysi)
 % hold on
 % step(dsysi)
 % return

 % feedback controller
 p1 = -100;
 p2 = -110;
 p3 = -200;
 p4 = -300;
 poles_k_c = [p1 p2 p3 p4]
 poles_k_d = exp(Ts.*poles_k_c)   
 % poles_k = exp(Ts*poles_k)   
 % return
 K_c = place(asysi.a,asysi.b,poles_k_c)
 K_d = place(dsysi.a,dsysi.b,poles_k_d)

 % Ki=K(4);                                
 % K=K(1:3)

 t = 0:0.001:.1;
 sys_cl_c = ss(Ai-bi*K_c,br,cTi,di)
 pole(sys_cl_c)
 step(sys_cl_c,t)
 hold on;
 % return

 sys_cl_d1 = c2d(sys_cl_c,Ts)
 step(sys_cl_d1,0:0.01:.1)
 pole(sys_cl_d1)

 sys_cl_d2 = ss(dsysi.a-dsysi.b*K_d,dsysi.b,dsysi.c,dsysi.d,Ts)
 pole(sys_cl_d2)
 step(sys_cl_d2,0:0.01:1)

First, you should use exp(Ts.*poles_k) rather than exp(Ts*poles_k)

Second, you need to use ss(.,Ts) for a discrete system rather than ss(.) which is for continuous systems.

As you may see in the result, the step responses of sys_cl_c and sys_cl_d1 are consistent. Still, the response of sys_cl_d2 doesn't match even though it has the same poles of sys_cl_d1

[Kevin]

so far i can follow your approaches but i guess my problem is based on the simulink model( v1 ). First of all i tried to realise the model by several blocks which obviously doesn't work for a discretized model, at least not as supposed. so maybe it's easier to take the state-space blocks( v2 ) instead but how to implement the integral and the feedback gain then? Again for the continuous case it's not the problem. the upper state space block represents the open loop system, the other state space block includes the observer gain L. firstly, L is added to bi, and then the entire observer is used for calculating the gains K and Ki.

 %system model
A = [0 1 0;
      0 -Bm/J k_phi/J;
      0 -k_phi/La -Ra/La];
b = [0 ; 0 ; 1/La];
cT = [1 0 0];
d = [0];
 asys=ss(A,b,cT,d)    ;  
 dsys=c2d(asys,Ts,'zoh') ;

 poles_obsv = [-1000 -1100 -2000];       
 L=place(asys.a',asys.c',poles_obsv);    
 L=L';

 % observer model
 Ah = A;
 bh = [b L];
 cTh = eye(3);
 dh = [0 0;0 0;0 0];
 asysh=ss(Ah,bh,cTh,dh);  

 %model for integral action
 Ai = [0 1 0 0;
    0 -Bm/J k_phi/J 0;
    0 -k_phi/La -Ra/La 0;
    1 0 0 0];
 bi = [0 ; 0 ; 1/La ; 0 ];
 bi = [bi [L; 0]];
 br = [0 0; 0 0 ; 0 0; -1 0];
 cTi = [1 0 0 0];
 di = [0];
 asysi=ss(Ai,bi,cTi,di);   

 % feedback controller
 poles_k_c = [-100 -110 -200 -300];
 K_c = place(asysi.a,asysi.b,poles_k_c)

 K = K_c(1,1:3)
 Ki = K_c(1,4)

what to do for K and Ki in the discretized case? - since you've already mentioned that gains differ between the discrete and continuous presentation. how to describe the observer together with the integral action then? - sry, but still i don't get it...

[Kevin]

ok, now i have a final proposal: look here

the upper state space block is the open loop model. below, you find an observer which is used for the calculation of a feedback gain and an integral action. based on a time continuous design the discretized version is derived. the related code would be:

 %system model
A = [0 1 0;
      0 -Bm/J k_phi/J;
      0 -k_phi/La -Ra/La];
b = [0 ; 0 ; 1/La];
cT = [1 0 0];
d = [0];
 asys=ss(A,b,cT,d)    ;  
 dsys=c2d(asys,Ts,'zoh') ;

 poles_obsv = [-1000 -1100 -2000];       
 L=place(asys.a',asys.c',poles_obsv);    
 L=L';

 % observer model
 Ah = A;
 bh = [b L];
 cTh = eye(3);
 dh = [0 0;0 0;0 0];
 asysh=ss(Ah,bh,cTh,dh);  
 dsysh=c2d(asysh,Ts,'zoh') ;

 %model for integral action
 Ai = [0 1 0 0;
    0 -Bm/J k_phi/J 0;
    0 -k_phi/La -Ra/La 0;
    1 0 0 0];
 bi = [0 ; 0 ; 1/La ; 0 ];
 bi = [bi [L; 0]];
 br = [0 0; 0 0 ; 0 0; -1 0];
 cTi = [1 0 0 0];
 di = [0];
 asysi=ss(Ai,bi,cTi,di);   
 dsysi=c2d(asysi,Ts,'zoh') ;

 % feedback controller
 poles_k_c = [-100 -110 -200 -300];
 K_c = place(asysi.a,asysi.b,poles_k_c)

 sys_cl = ss(Ai-bi*K_c,br,cTi,di);
 sys_cl_d1 = c2d(sys_cl,Ts,'zoh');

would you agree on that solution? at least the result looks pretty sufficient...

[Babak]

Looks fine... Again, remember that exp(Ts.*poles_k) and exp(Ts*poles_k) are different! You need to use exp(Ts.*poles_k) to convert your continuous system poles to the discrete. P.S.: Due to my limited internet access here I cannot view the links you've attached.

[Azzi Abdelmalek]

Babak what is the difference between 2*[1 2 3] and 2.*[1 2 3]? the result is the same

Ts is the sample time

Answer 4

[Azzi Abdelmalek]

% Kevin, let us compare for stable system. we will obtain the same result
close all
Ts=0.001

% stable continuous system

 N=1;D=poly([-12 -12 -12])
[A,b,cT,d]=tf2ss(N,D)
 asys=ss(A,b,cT,d)    ;

% initial system + integral action (there is an error in your model)

Ai = [[A   zeros(3,1)];-cT  0]   %-cT instead cT
bi = [b; 0 ];
br = [zeros(3,1);1];   % 1 instead -1
cTi = [cT 0];
di = [0];
asysi=ss( Ai,bi,cTi,di);

% controller parameters

poles_k =[-100  -110  -200  -300];
K = place(asysi.a,asysi.b,poles_k);    
Ki=K(4);K=K(1:3);

% simulation

 t = 0:Ts:.1;
sys_cl = ss(Ai-bi*[K Ki],br,cTi,di);
step(sys_cl,t)

%discrete system . here was the trick

sys_cld=c2d(sys_cl,Ts);  % 'zoh' is the default setting

%you cant use

%sys_cld = ss(dsysi.a-dsysi.b*[Kd Kdi],br,dsysi.c,dsysi.d);

%discret system simulation

figure
step(sys_cld,t)

[Azzi Abdelmalek]

to apply to your motor

close all
Ra  = 2.06;            
La = 0.000238;         
k_phi = 1/((406*2*pi)/60); 
J = 1.07e-6;           
Bm = 12e-7;   

Ts=0.002;

%general model %states: theta, dtheta, i, derr

A = [0 1 0;
  0 -Bm/J k_phi/J;
  0 -k_phi/La -Ra/La];
b = [0 ; 0 ; 1/La];
cT = [1 0 0];
d = [0];
asys=ss(A,b,cT,d);           
dsys=c2d(asys,Ts,'zoh');

%model for integral action

Ai = [[A   zeros(3,1)];-cT  0]
bi = [b; 0 ];
br = [zeros(3,1); 1];
cTi = [cT 0];
di = [0];
asysi=ss(Ai,bi,cTi,di);

% feedback controller

p1 = -100;
p2 = -110;
p3 = -200;
p4 = -300;
poles_k = [p1 p2 p3 p4];                
K = place(asysi.a,asysi.b,poles_k);    
Ki=K(4);                                
K=K(1:3);                               

sys_cl = ss(Ai-bi*[K Ki],br,cTi,di);
t = 0:Ts:0.1;
step(sys_cl,t)

%discret system hold on

sys_cld = c2d(sys_cl,Ts)
step(sys_cld,t)