Good and valued: Solution ver. 1,
clc
clear all
close all
%% consider a 2 input, 2 output system,
nstates = 4;
ninputs = 2;
noutputs = 2;
A1=[0 0 -178 178;0 0 590 -500;-256 256 0 0; 0 -1 0 0.0000666];
B1=[8 -5;48 -63;-0.01 2.5; 0 0];
C1=[1 0 0 0;0 0 0 1];
D1=[0 0; 0 0];
Ts = 10e-3; % Sampling time in s
% Aa=[0 0 -178 178;0 0 590 -500;-256 256 0 0; 0 -1 0 0.0000666];
% Ba=[8 -5;48 -63;-0.01 2.5; 0 0];
% Ca=[1 0 0 0;0 0 0 1];
% Da = [0 0; 0 0];
% % you can use 'rss' and 'ssdata' to get random A, B, C, D matrices
% sys = rss(2,2,2) % to generate a random 2x2 system
% [A,B,C,D] = ssdata(sys); % to provide A, B, C, D matrices for this system
[num1 den1] = ss2tf(A1,B1,C1,D1,1) % iu = 1
[num2 den2] = ss2tf(A1,B1,C1,D1,2); % iu = 2
%First input to first output only:
sys1 = tf(num1(1,:),den1) %contin
%Second input to second output only:
sys2 = tf(num2(2,:),den2); %contin
%First input to second output only:
sys3 = tf(num1(2,:),den1); %continuous
%Second input to first output only:
sys4 = tf(num2(1,:),den2); %continuous
[A, B, C, D] = tf2ss(num1(1,:),den1)
asys = ss(A,B,C,D)
opt2 = c2dOptions('Method','tustin','ThiranOrder',4,'DelayModeling','state');
%dsys=c2d(asys,Ts,opt2)
% observer model
%poles_obsv = exp(Ts*[-700 -310 -100 -90]);
poles_obsv = [-700 -310 -100 -90];
L=place(asys.a',asys.c',poles_obsv);
L=L';
Ah = asys.A;
bh = [asys.B L];
%cTh = eye(nstates);
%dh = [0 0;0 0;0 0;0 0];
%asysh=ss(Ah,bh,cTh,dh);
%figure
%step(asysh)
%add a state variable for the integral of the output error. This has the effect of adding an integral term to our controller which is known to reduce steady-state error.
%model for integral action
Ai = [[asys.A zeros(nstates,1)];-asys.C 0];
bi = [asys.B;0];
br = [zeros(nstates,1);1];
ci = [asys.C 0];
di = [0];
asysi=ss(Ai,bi,ci,di);
%Other augmentation scheme:
% Plant augmentation
Aaug=[asys.A zeros(nstates,1); zeros(1,nstates+1)];
Aaug(nstates+1,nstates)=1;
Baug=[asys.B;0];
Caug=eye(nstates+1);
Daug=zeros(nstates+1,1);
Plantcs=ss(Aaug,Baug,Caug,Daug);
% feedback controller
p1 = -800 + 800i;
p2 = -800 - 800i;
p3 = -400 - 400i;
%p1 = -110;
%p2 = -310;
%p3 = -500;
p4 = -400 + 400i;
p5 = -90;
%poles_k = exp(Ts*[p1 p2 p3 p4 p5]);
poles_k = [p1 p2 p3 p4 p5];
K = place(asysi.a,asysi.b,poles_k)
Ko = place(asys.A,asys.B,[p1 p2 p3 p4])
s4 = size(asys.A,1);
Z = [zeros([1,s4]) 1];
N = (1\([asys.A,asys.B;asys.C,asys.D]))*Z';
Nx = N(1:s4);
Nu = N(s4+1);
Nbar=Nu + Ko*Nx
%observer design:
At = [ asys.A-asys.B*Ko asys.B*Ko
zeros(size(asys.A)) asys.A-L*asys.C ];
Bt = [ asys.B*Nbar
zeros(size(asys.B)) ];
Ct = [ asys.C zeros(size(asys.C)) ];
%If you'd like to eliminate steady state error as much use Nbar:
% compute Nbar
%poles_k_d = exp(Ts.*poles_k)
%K_d = place(dsysi.a,dsysi.b,poles_k_d)
Ki=K(nstates+1);
K=K(1:nstates);
s4 = size(asys.A,1);
Z = [zeros([1,s4]) 1];
N = (1\([asys.A,asys.B;asys.C,asys.D]))*Z';
Nx = N(1:s4);
Nu = N(s4+1);
Nbarq=Nu + K*Nx
%Well performing (adjusted) system:
syso = ss(At,Bt,Ct,0);
isstable(syso)
syso2 = minreal(syso)
%isminphase(syso)
f = 3;
t = 0:Ts:f;
figure
step(syso2,[0 f])
title('Continuous With observer flat(cold) start')
sysod = c2d(syso2,Ts,opt2)
isstable(sysod)
u = 0.001*ones(size(t));
x0 = [0.01 0 0 0];
figure
lsim(sysod,zeros(size(t)),t,[x0]);
title('Discrete sampled with observer and conditions hot start')
xlabel('Time (sec)')
ylabel('Output')
res_agent = sysod;
figure
w = logspace(2,6,1000);
bode(res_agent,w),grid
figure
rlocus(res_agent)
figure
nyquist(res_agent)
ms = minreal(res_agent);
Gcc = ms;
zms = zero(ms) % closed-loop zeros
pms = pole(ms)
tsim = 1;
setpointamp = 1; %
setpointapptime = 0.001; %
defaultsetpointpos = 0; %
Conf = RespConfig(Bias=defaultsetpointpos,Amplitude=setpointamp,Delay=setpointapptime)
%figure
%step(Gcc,[0 tsim], Conf)
%title('Plant+Controller Closed Loop system step response')
[wngcc,zetagcc,pgcc] = damp(Gcc)
sigma(Gcc)
margin(Gcc)
diskmargin(Gcc)
stepinfo(Gcc)
disp('Whether the system is stable, minimum phase and proper:')
isstable(sysod)
isminphase(tfdata(tf(sysod),'v'))
isproper(sysod)
Qc= ctrb(A,B)
controllab = rank(Qc)
hasdelay(sysod)
%tdc3 = totaldelay(sysod)
%sysndc3 = absorbDelay(sysod)
isallpass(tfdata(tf(sysod),'v'))
%Less important:
%sys_cl = ss(Ai-bi*[K Ki],br,ci,di)
%[a,b] = ss2tf(sys_cl.A,sys_cl.B,sys_cl.C,sys_cl.D)
%[A,B,C,D] = tf2ss(Nbar*a,b)
%sys_cl = ss(A,B,C,D)
%figure
%step(sys_cl,t)
%discrete system
%hold on
%sys_cld = c2d(sys_cl,Ts)
%figure
%step(sys_cld)
%title('Discrete sampled System')
%In similar approach to continue for 1-2, 2-1, 2-2 SISO parts of the MIMO.
%Contact me more if you need more tailored, specific service, solution,
%product development.
%Constructed by
%https://independent.academia.edu/PMazniker
%+380990535261, https://join.skype.com/invite/oXnJhbgys7oW
%https://diag.net/u/u6r3ondjie0w0l8138bafm095b
%https://github.com/goodengineer
%https://orcid.org/0000-0001-8184-8166
%https://willwork781147312.wordpress.com/portfolio/cp/
%https://www.youtube.com/channel/UCC__7jMOAHak0MVkUFtmO-w
%https://nanohub.org/members/130066
%https://pangian.com/user/hiretoserve/
%https://substack.com/profile/191772642-paul-m
