Digital Control Octave Scripts - Revision history

Frohro: /* Kalman Filter RC Circuit Example */

2016-05-26T21:38:21Z

Kalman Filter RC Circuit Example

← Older revision		Revision as of 14:38, 26 May 2016
Line 61:		Line 61:
	C = 1000e-6;		C = 1000e-6;
	A = -1/(R*C);		A = -1/(R*C);
	B = 1000;%/(R*C);		B = 1000/(R*C);
	C_cont = 1;		C_cont = 1;
	cont_sys = ss(A,B,C_cont);		cont_sys = ss(A,B,C_cont);
Line 74:		Line 74:
	[y,t,x] = lsim(disc_sys,u) ; % The original discretized system without any noise.		[y,t,x] = lsim(disc_sys,u) ; % The original discretized system without any noise.
	w = .01*randn(size(t)); % random state noise		w = .01*randn(size(t)); % random state noise
	v = sqrt(.01)*randn(size(t)); % random measurement noise		v = sqrt(.1)*randn(size(t)); % random measurement noise
	% Note: to add w to x and retain u, u = u + 1/Gamma*w		% Note: to add w to x and retain u, u = u + 1/Gamma*w
	[yn,t,xn] = lsim(disc_sys,u+1/Gamma*w); %noisy due to state errors.		[yn,t,xn] = lsim(disc_sys,u+1/Gamma*w); %noisy due to state errors.
Line 81:		Line 81:
	R = cov(v,v')		R = cov(v,v')

	% ~~First: show how to use~~ the built in ~~kalman function to get steady state results.~~		% Use the built in dlqe (Discrete Linear Quadratic Kalman Estimator) fro comparison.

	% Block Diagram
	~~% u +-------+ ^~~
	~~% +---------------------------->\| \|-------> y~~
	~~% \| +-------+ + y \| est \| ^~~
	~~% u ----+--->\| \|----->~~(+)~~------>\| \|-------> x~~
	~~% \| sys \| ^ + +-------+~~
	~~% w -------->\| \| \|~~
	~~% +-------+ \| v~~
	%
	~~% Q = cov (w, w') R = cov (v, v') S = cov (w, v')~~

	[m,p,z,e] = dlqe(Phi,[],disc_sys.c,Q,R)		[m,p,z,e] = dlqe(Phi,[],disc_sys.c,Q,R)

Frohro: /* Kalman Filter RC Circuit Example */

2016-05-26T18:51:57Z

Kalman Filter RC Circuit Example

← Older revision		Revision as of 11:51, 26 May 2016
Line 52:		Line 52:
	This is for a series RC circuit. It shows how to do it yourself using equations we derived in class and using dlqe().		This is for a series RC circuit. It shows how to do it yourself using equations we derived in class and using dlqe().
	<nowiki>		<nowiki>


	% This is a demo of the Kalman filter on an RC circuit.		% This is a demo of the Kalman filter on an RC circuit.
	% Rob Frohne, for ENGR 454, Digital Control Systems.		% Rob Frohne, for ENGR 454, Digital Control Systems.

Frohro: /* System Identification Simulation */

2016-05-26T18:51:29Z

System Identification Simulation

@@ Line 46: / Line 46: @@
 plot(t,y,t,y_fit)
 sum_err = (y-y_fit)'*(y-y_fit)
 </nowiki>

Frohro: Created page with "===System Identification Simulation=== % Digital Control System Identification Example pkg load control clear all; close all; % Set up the system to fit to. A=diag([..."

2016-04-01T01:09:09Z

Created page with "===System Identification Simulation=== <nowiki> % Digital Control System Identification Example pkg load control clear all; close all; % Set up the system to fit to. A=diag([..."

New page

===System Identification Simulation===
<nowiki>
% Digital Control System Identification Example
pkg load control
clear all;
close all;
% Set up the system to fit to.
A=diag([-1,-2]) % Make sure the system is stable.
B=[1,1]'
C=[1,2]
%A=diag(-1)
%B=1
%C=1
sys = ss(A,B,C); % sys is the system to fit.
% Move to the sampled domain.
T=0.01; % T is the sampling interval.
N=10000; % N is the number of training points.
t=0:T:N*T;
sys_discrete = c2d(sys,T);
sys_discrete_tf = tf(sys_discrete)
if(!isstable(sys))
error('System is not stable!\n') % Check for stability.
end
if(!isctrb(sys)) % Check to make sure the system is contrololable,
error('System is not controllable!\n')
end
if(!isobsv(sys)) % and observable. If it is not you need to fix it!
error('System is not observable!\n')
end
% Now set up the system identification.
u=randn(size(t));
[y,t]=lsim(sys,u,t);
% The order of the 0fitted system is n.
n=2
Y=[[y(n+1:N+1)]]; % Note: Octave's lowest index is 1, not 0.
for k=n:N
PSI(k-n+1,:)=[y(k:-1:k-n+1)', u(k:-1:k-n+1)];
end
PSI;
THETA=(PSI'*PSI)^(-1)*PSI'*Y
sys_fit=tf(THETA(n+1:2*n),[1,-THETA(1:n)'],T) % tf(num,den,T)
[y_fit,t]=lsim(sys_fit,u,t);
err=(y-y_fit);
plot(t,err)
figure;
plot(t,y,t,y_fit)
sum_err = (y-y_fit)'*(y-y_fit)
</nowiki>