MàJ de script MATLAB
This commit is contained in:
parent
3058e66faf
commit
eefb551d9c
1 changed files with 142 additions and 0 deletions
142
identificationScript.m
Executable file
142
identificationScript.m
Executable file
|
|
@ -0,0 +1,142 @@
|
|||
%% Initialisation
|
||||
clc
|
||||
clear all
|
||||
close all
|
||||
%load expe1.mat simout
|
||||
%load expe2_multisine.mat
|
||||
%load expe3_multisine.mat
|
||||
%load expe4_multisine.mat
|
||||
load expe5_multisine.mat
|
||||
|
||||
%% Variables
|
||||
N = 2^9; % Nombre
|
||||
Te = 5e-2; % Temps d'échantillonage
|
||||
fdeb = 0; % fréquence de début d'échantillonage
|
||||
ffin = 4; % fréquence de fin d'éch.
|
||||
BP = [fdeb, ffin]; % Bande passante
|
||||
REV = false; % Si on passe de fdeb à ffin ou l'inverse (booléen)
|
||||
SHOW = true; % Montre les résultats (graphes) du MLBS (booléen)
|
||||
|
||||
%[u, t] = insapack.mlbs(N, Te, BP, REV, SHOW);
|
||||
[u, t] = insapack.multisine(N, Te, BP, false,'all',REV, SHOW);
|
||||
u = 4*u;
|
||||
% var.time = t.';
|
||||
% var.signals.values=u.';
|
||||
% var.signals.dimensions= 1 ;
|
||||
%%
|
||||
% Muligens litt for hoey frekvens, den maa kanskje skrus litt ned. Deretter
|
||||
% finner vi modell ved bruk av div. verktoey.
|
||||
|
||||
% figure()
|
||||
% hold on;
|
||||
% plot(simout.time, simout.signals.values)
|
||||
% plot(var.time, var.signals.values)
|
||||
|
||||
%%% N4SID
|
||||
tn = simout.time; % On enlève les dix primières valeurs
|
||||
un = var.signals.values; % On enlève les dix primières valeurs
|
||||
yn = simout.signals.values; % On enlève les dix primières valeurs
|
||||
data = iddata(yn,un,Te);
|
||||
ordre = 2;
|
||||
sys = n4sid(data, ordre, 'Ts',Te);
|
||||
y = lsim(ss(sys), un, tn);
|
||||
H = tf(sys);
|
||||
G = (H)/(1-H);
|
||||
numG = G.Numerator;
|
||||
denumG = G.Denominator;
|
||||
|
||||
|
||||
%%% Freq
|
||||
[f0,U0,Y0,G0,sigU2,sigY2,sigUY2,sigG2,b,rho] = insapack.non_param_freq(un,yn,Te);
|
||||
w0 = 2*pi*f0;
|
||||
|
||||
%%% Identification via N4SID in frequency-domain
|
||||
data = iddata(Y0,U0,Te,'Frequency',w0); %Ici peut-etre changer avec f0, comme ça il marche très mieux avec Loewner
|
||||
Hn4sid_fd = n4sid(data,ordre);
|
||||
y2 = lsim(ss(Hn4sid_fd), un, tn);
|
||||
|
||||
%%% Identification via Loewner
|
||||
wid = 2*pi*f0(f0<BP(end)/2);
|
||||
wRange = 1:floor(length(wid)/2)*2;
|
||||
wid = 2*pi*f0(wRange);
|
||||
G0id = G0(wRange);
|
||||
[la,mu,W,V,R,L] = insapack.data2loewner(wid,G0);
|
||||
opt.target = ordre;
|
||||
[hr,info] = insapack.loewner_tng(la,mu,W,V,R,L,opt);
|
||||
Hloe = dss(info.Ar,info.Br,info.Cr,info.Dr,info.Er);
|
||||
Hloe = stabsep(Hloe);
|
||||
y3 = lsim(ss(Hloe), un, tn);
|
||||
|
||||
%%% Plots
|
||||
figure()
|
||||
subplot(211)
|
||||
plot(tn, un)
|
||||
legend("Signal d'entrée (d'exitation)")
|
||||
subplot(212)
|
||||
hold on;
|
||||
plot(tn, yn)
|
||||
plot(tn, y)
|
||||
plot(tn, y2)
|
||||
plot(tn, y3)
|
||||
legend('Signal de sortie observée','Sortie simuléee avec n4sid en temporel', 'Sortie simuléee avec n4sid en fréquenciel', 'Sortie simuléee avec Loewner')
|
||||
|
||||
|
||||
var2 = var;
|
||||
var2.signals.values(1:100) = 1;
|
||||
var2.signals.values(101:200) = 3.3;
|
||||
var2.signals.values(201:300) = 0;
|
||||
var2.signals.values(301:400) = 1.5;
|
||||
var2.signals.values(401:503) = 2.7;
|
||||
|
||||
%% Comment est-ce qu'on a trouvé la valeur du P dans notre système rail
|
||||
Gc=d2c(G)
|
||||
consigne = 1;
|
||||
P = [ 0.8 1 1.5 2 3 5];
|
||||
figure()
|
||||
hold on;
|
||||
|
||||
for i = 1:1:length(P)
|
||||
Gcf = P(i)*Gc/(1+P(i)*Gc);
|
||||
step(Gcf,10)
|
||||
end;
|
||||
yline(0);
|
||||
title("Reponse d'échelon du système rail");
|
||||
legend("Kp=0.8","Kp=1","Kp=1.5","Kp=2","Kp=3","Kp=5");
|
||||
|
||||
|
||||
%% Bille
|
||||
g = 9.81
|
||||
|
||||
a = 10
|
||||
b = 0.2 % V/cm sur le rail, +-10V pour le rail de +-50cm
|
||||
Kb = b*g
|
||||
T = 0.22
|
||||
Kp = 1
|
||||
|
||||
Gb = tf([Kb], [1 0 0])
|
||||
|
||||
figure()
|
||||
bode(Gb)
|
||||
title('Bode système rail')
|
||||
|
||||
Cb = tf([a*T 1], [T, 1])*Kp; %Correcteur d'avance de phase
|
||||
|
||||
figure()
|
||||
bode(Cb)
|
||||
title('Bode correcteur avance de phase')
|
||||
|
||||
%% marge de phase
|
||||
figure()
|
||||
title("Marge de phase et gain pour le système complet");
|
||||
Hc = d2c(H);
|
||||
Lbo = Cb*Hc*Gb;
|
||||
Lbf = Lbo / (1+Lbo);
|
||||
Margin = allmargin(Lbf);
|
||||
h=bodeplot(Lbf);
|
||||
margin(Lbf);
|
||||
title({'Marge de phase et gain pour le système complet','Gm= Inf Pm =-40.6 deg (at 5.97rad/s)'});
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
Loading…
Reference in a new issue