State space model with disturbance

@Paul If I have a system x_dot= [500 0;0 0.1]x+[50;0]u+[0 0.001]d and y= [0 1]x. So I want to make an LQR controller but my problem with the disturbance term and how can I deal with it?

[Sam Chak]

@Mohammed Yakoob

Do you have the knowledge about the disturbance d?

Is the disturbance d bounded?

Try providing these info before @Paul comes.

[Mohammed Yakoob]

Thanks for your comments @Sam Chak!! I think The disturbance will be bounded because we it will be the outputs of the electric circuit with different loads for example first type will be non linear and the second will be harmonic and so on!!

Answer 1

[Sam Chak]

Hi @Mohammed Yakoob

I refer to the state-space system from another post that you opened in the link:

https://www.mathworks.com/matlabcentral/answers/1684219-non-linear-loads

where u is the controller, and w is the unmatched bounded disturbance that you mentioned above.

If the unmatched disturbance w is measuarable, then the proposed controller u is given by

where ζ and Ω are the two control parameters to be determined. Note that ζ is damping ratio and Ω is the undamped natural frequency of the control system.

If you let the autotuner to design the Linear Quadratic Regulator (LQR) for you and obtain the gains , then you can apply

.

Try it and show the result to see if it works or not.

[Mohammed Yakoob]

Hi Sam Chak and thanks for your reply!! My problem is how can I get the disturbance from these different loads to return them as disturbance. I tried something in the following code but I don't know if I correct or not. Also, if you ask me why don't take the load from electrical circuit that because I don't know how to do that and connect the loads to my model

[Mohammed Yakoob]

My codes : clc close all clear %% Plant Paramiters f= 60; Ts= 1/f; %% one period Th=Ts/2; %%half period Wo= 2*pi*f; Vdc= 300; C_st=15*10^-6; L_t=2*10^-3; R_l= 40; %resistance load R_li=3; %resistance of the line

%% Implement a Single Phase Microgrid in a contious time form % x_dot= AX+BU+dW ; Y= CX+DU+0W. % X=[il Vg]' ; U=[Vsw] ; d= [ig]; Y=[Vg]. % for i= 1:1000 Aop=[0 1/L_t; 1/C_st 0]; Bop= [1/L_t; 0]; d=[0; -1/C_st]; Cop=[0 1]; Dop= [0]; B= [Bop d]; %% Imlement the system sys=ss(Aop,B, Cop, Dop); sys_d= c2d(sys,0.0001); [Ad,Bd,Cd,Dd]=ssdata(sys_d); %% Performance of The Open Loop System t = 0.000:0.0001:Th; V_r =Vdc*sin(Wo* t); %% Cotrollability of the system Co= ctrb(sys); rank(Co); %% LQR_LMI Ac=Aop; Bc= Bop; C= Cop; nx = size(Ac,2); nu = size(Bc,2); % closed loop % LQR weights Q = diag([1 1e-5]); R = 1; Klqr = lqr(Ac,Bop,Q,R,[]);

% % % LMI % % A=Ac; % % B1 = eye(nx); % % B2=Bc; % % % % C1 = eye(nx); % % D11 = zeros(nx,nx); % % D12 = zeros(nx,nu); % % % % C2 = [sqrt(Q); % % zeros(nu,nx)]; % % D21 = zeros(nx+nu,nx); % % D22 = [zeros(nx,nu); % % sqrt®]; % % % % Pp = ltisys(A, [B1 B2], [C1; C2], [D11 D12; D21 D22]); % % % % siz = size(D22); % % gamma0 = 0; % % obj = [gamma0 0 0 1]; % % % % [gopt,h2opt,Klmi,Pcl,X] = msfsyn(Pp,siz,obj); % % % % Klqr % % Klmi %% Closed loop system Acl= [Aop-(Bop*Klqr)]; Bcl= [B]; Ccl= [0 1]; Dcl= [0 0]; states = {'i_L' 'V_g'}; inputs1 = {'V_r' 'dis'}; outputs1 = {'V_g'}; poles= eig(Acl); sys_cl = ss(Acl,Bcl,Ccl,Dcl,'statename',states,'inputname',inputs1,'outputname',outputs1); %% At the disturebance eqaul to zero input1 = [V_r; zeros(size(t))]; %%assume tha the disturebance =0. [y1,t,x]=lsim(sys_cl,input1,t); figure() plot(t,y1) ylabel('V_out') xlabel('time') title('LQR controller Without the Disterbance Term') grid on %% LQR when the system has a load term at R=40 Ohm dis = (V_r/R_l+R_li); input2 = [V_r' dis'.*0.01]; [y2,t,x]=lsim(sys_cl,input2,t); figure() plot(t,y2) ylabel('V_out') xlabel('time') title('LQR controller With the Disterbance Term at R=40Ohm') grid on %% 1) PERFORMANCE AGAINST UNKNOWN LOAD z1=complex(0,-Wo*62.86e-6); z2= complex(0.35,Wo*223.8e-3); z3= complex(228); z4=complex(76,-Wo*10e-6); z5=complex(152,Wo*111.9e-3); Z11= (1/z1)+ (1/z2) +(1/z3); Z22=(1/z4)+(1/z5);

n= numel(t);
% %  for  m= 1:n 
% %     
% %  if m<= 36
% %          dis11= (y1(m,1)/R_li+Z11)';
% %          dis11(1,m)=dis11;
% %  elseif m<=n-36
% %      
% %         dis12= (y1(36+m,1)/R_li+Z11+Z22)';
% %         dis12(1,m)=dis12;
% %     %  end
% %      
% %  end 
%      dis1=[dis11 dis12];
for  m= 1:n 

 if m<= 36
         dis11= 1e-1*[cos(5*2*pi*m')+cos(1*2*pi*m')+cos(.1*2*pi*m')];
         dis11(1,m)=dis11;
 elseif m<=n-36

        dis12= 1*[cos(4*2*pi*m')+cos(2*2*pi*m')+cos(.2*2*pi*m')];
        dis12(1,m)=dis12;

 end

 end 
dis1=[dis11 dis12];
 input3 = [V_r' dis1'];
 [y3,t,x]=lsim(sys_cl,input3,t);
 figure()
 plot(t,y3)
 ylabel('V_out')
 xlabel('time')
 title('PERFORMANCE AGAINST UNKNOWN LOAD')
 grid on

 %%(2) PERFORMANCE AGAINST HARMONIC LOAD
 dis = ((V_r/R_l+R_li)+7*sin(Wo*t'));
input3 = [V_r' dis'.*0.001];
[y3,t,x]=lsim(sys_cl,input3,t);
figure()
plot(t,y3)
ylabel('V_out')
xlabel('time')
title('PERFORMANCE AGAINST HARMONIC LOAD')
grid on

[Sam Chak]

@Mohammed Yakoob Can you edit your post so that your MATLAB code is in readable format?

[Sam Chak]

@Mohammed Yakoob

Can you confirm if the following LQR u?

Can you confirm which of the following disturbance is true?

... as displayed in your code

... according to Ohm's Law

where .

The lsim() function can take multiple inputs, but the first input and the second input are fed into different "systems". You can find the desired output using the superposition principle if it holds.

Else, can you consider using the more direct approach with the ode45() method to feed in the LQR u and the disturbance w?

with initial condition , and .

Answer 2

[Mohammed Yakoob]

clc

close all

clear

%% Plant Paramiters

f= 60;

Ts= 1/f; %% one period

Th=Ts/2; %%half period

Wo= 2*pi*f;

Vdc= 300;

C_st=15*10^-6;

L_t=2*10^-3;

R_l= 40; %resistance load

R_li=3; %resistance of the line

%% Implement a Single Phase Microgrid in a contious time form

% x_dot= AX+BU+dW ; Y= CX+DU+0W.

% X=[il Vg]' ; U=[Vsw] ; d= [ig]; Y=[Vg].

% for i= 1:1000

Aop=[0 1/L_t; 1/C_st 0];

Bop= [1/L_t; 0];

d=[0; -1/C_st];

Cop=[0 1];

Dop= [0];

B= [Bop d];

%% Imlement the system

sys=ss(Aop,B, Cop, Dop);

sys_d= c2d(sys,0.0001);

[Ad,Bd,Cd,Dd]=ssdata(sys_d);

%% Performance of The Open Loop System

t = 0.000:0.0001:Th;

V_r =Vdc*sin(Wo* t);

%% Cotrollability of the system

Co= ctrb(sys);

rank(Co);

%% LQR_LMI

Ac=Aop;

Bc= Bop;

C= Cop;

nx = size(Ac,2);

nu = size(Bc,2);

% closed loop

% LQR weights

Q = diag([1 1e-5]);

R = 1;

Klqr = lqr(Ac,Bop,Q,R,[]);

% % % LMI

% % A=Ac;

% % B1 = eye(nx);

% % B2=Bc;

% %

% % C1 = eye(nx);

% % D11 = zeros(nx,nx);

% % D12 = zeros(nx,nu);

% %

% % C2 = [sqrt(Q);

% % zeros(nu,nx)];

% % D21 = zeros(nx+nu,nx);

% % D22 = [zeros(nx,nu);

% % sqrt(R)];

% %

% % Pp = ltisys(A, [B1 B2], [C1; C2], [D11 D12; D21 D22]);

% %

% % siz = size(D22);

% % gamma0 = 0;

% % obj = [gamma0 0 0 1];

% %

% % [gopt,h2opt,Klmi,Pcl,X] = msfsyn(Pp,siz,obj);

% %

% % Klqr

% % Klmi

%% Closed loop system

Acl= [Aop-(Bop*Klqr)];

Bcl= [B];

Ccl= [0 1];

Dcl= [0 0];

states = {'i_L' 'V_g'};

inputs1 = {'V_r' 'dis'};

outputs1 = {'V_g'};

poles= eig(Acl);

sys_cl = ss(Acl,Bcl,Ccl,Dcl,'statename',states,'inputname',inputs1,'outputname',outputs1);

%% At the disturebance eqaul to zero

input1 = [V_r; zeros(size(t))]; %%assume tha the disturebance =0.

[y1,t,x]=lsim(sys_cl,input1,t);

figure()

plot(t,y1)

ylabel('V_out')

xlabel('time')

title('LQR controller Without the Disterbance Term')

grid on

%% LQR when the system has a load term at R=40 Ohm

dis = (V_r/R_l+R_li);

input2 = [V_r' dis'.*0.01];

[y2,t,x]=lsim(sys_cl,input2,t);

figure()

plot(t,y2)

ylabel('V_out')

xlabel('time')

title('LQR controller With the Disterbance Term at R=40Ohm')

grid on

%% 1) PERFORMANCE AGAINST UNKNOWN LOAD

z1=complex(0,-Wo*62.86e-6);

z2= complex(0.35,Wo*223.8e-3);

z3= complex(228);

z4=complex(76,-Wo*10e-6);

z5=complex(152,Wo*111.9e-3);

Z11= (1/z1)+ (1/z2) +(1/z3);

Z22=(1/z4)+(1/z5);

n= numel(t);

% % for m= 1:n

% %

% % if m<= 36

% % dis11= (y1(m,1)/R_li+Z11)';

% % dis11(1,m)=dis11;

% % elseif m<=n-36

% %

% % dis12= (y1(36+m,1)/R_li+Z11+Z22)';

% % dis12(1,m)=dis12;

% %

% % end

% %

% % end

% dis1=[dis11 dis12];

for m= 1:n

if m<= 36

dis11= 1e-1*[cos(5*2*pi*m')+cos(1*2*pi*m')+cos(.1*2*pi*m')];

dis11(1,m)=dis11;

elseif m<=n-36

dis12= 1*[cos(4*2*pi*m')+cos(2*2*pi*m')+cos(.2*2*pi*m')];

dis12(1,m)=dis12;

end

end

dis1=[dis11 dis12];

input3 = [V_r' dis1'];

[y3,t,x]=lsim(sys_cl,input3,t);

figure()

plot(t,y3)

ylabel('V_out')

xlabel('time')

title('PERFORMANCE AGAINST UNKNOWN LOAD')

grid on

%% (2) PERFORMANCE AGAINST HARMONIC LOAD

dis = ((V_r/R_l+R_li)+7*sin(Wo*t'));

input3 = [V_r' dis'.*0.001];

[y3,t,x]=lsim(sys_cl,input3,t);

figure()

plot(t,y3)

ylabel('V_out')

xlabel('time')

title('PERFORMANCE AGAINST HARMONIC LOAD')

grid on