Het Tauforum

Code: Select all

close all
clear all
clc
% Het systeem is simpel:
%    twee massa's (parallel) verbonden met een veer en een demper.
%    beide massa's hebben wrijving met de grond
%    de eerste massa wordt aangedreven middels een kracht.
% Het probleem van dit voorbeeld zou een positioneringsprobleem kunnen
% zijn, waarbij:
%    de tweede massa moet gepositioneerd worden
%    de acceleratie van de eerste/tweede massa binnen bepaalde grenzen moet
%    blijven, maar bijv. ook de veer+demper kracht
%    (en deze moeten dus als output moeten worden opgenomen)
%
% Tekenafspraak: indrukken van de veer/damper levert een POSITIEVE kracht

% definieer eerst alle te gebruiken symbolen
syms ddz1 dz1 z1 ddz2 dz2 z2      % alle coordinaten (+afgeleiden)
syms u1                           % alle inputs
syms m1 m2 k d d1 d2              % alle paramters

% definieer tussenvergelijkingen (alle krachten)
F_veer   =  (z1 -z2 ) * k;
F_damper =  (dz1-dz2) * d;
F1_wrijf =   dz1      * d1;
F2_wrijf =   dz2      * d2;

% stel de bewegingsvergelijkingen op
ddz1 = (u1 - F_veer - F_damper - F1_wrijf)/m1;
ddz2 = (     F_veer + F_damper - F2_wrijf)/m2;

% definieer een toestand en input/output vector voor state space
x = [dz1 ; z1 ; dz2 ; z2];
u = [u1];

% stel nu f(x,u) samen
fx = [ddz1 ; dz1 ; ddz2    ; dz2     ];     % volgens: x' = f(x,u)
hx = [  z1 ;  z2 ; m1*ddz1 ; m2*ddz2 ; F_veer+F_damper];                           % volgens: y  = h(x,u)

% dan nu de A,B,C,D matrix berekenen
A = jacobian(fx,x);
B = jacobian(fx,u);
C = jacobian(hx,x);
D = jacobian(hx,u);

% definieer de waarden voor de parameters!
k  = 10;
d  =  1;
d1 =  5;
d2 = 10;
m1 = 10;
m2 =  5;

% vul deze waarden in. en converteer deze van symbolic > double
A = double(eval(A));
B = double(eval(B));
C = double(eval(C));
D = double(eval(D));

% maak er een state-space van
sys = minreal(ss(A,B,C,D));

% en plot de afzonderlijke bode-diagrammen (input,output)=(1,2) betekent:
%    laat de positie van massa 2 zien onder invloed van input 1, de kracht
input_output  = [
    1 1 ;   % welke input wil je zien?  [ 1 .. length(u)  ]
    1 2 ;   % welke output wil je zien? [ 1 .. length(hx) ]
];
for i=1:length(input_output)
    bode(sys(input_output(2,i),input_output(1,i))); hold on;
end

% we kunnen nu een stap op de ingang zetten van 50 N (tussen .1 en 2 s.):
dt = 1e-3;
t  = [0:5000]*dt;
ut = zeros(size(t)); ut( find( (.1<t) & (t<2) ) )=50;
yt = lsim(sys,ut,t);

figure
subplot(311); plot(t,[yt(:,1) yt(:,2) ut(:)./10]); legend('z_1','z_2','0.1*F(t)',4);
subplot(312); plot(t,yt(:,2)-yt(:,1));             legend('z_2-z_1',0);
subplot(325); plot(t,[yt(:,3) yt(:,4)]);           legend('z"_1','z"_2',0);
subplot(326); plot(t,[yt(:,5)]);                   legend('F_{veer}+F_{demper}',0);