1

Simple oscillator

(oscillator.mdl) Here we obtain a Simulink model of the simple spring-mass-damper system described by my + cy + ky = 0 To obtain a Simulink model, we rst note that y = k c y y m m

Simple oscillator (oscillator)

yddot

1 s

ydot

1 s

y Scope

0.1 c/m 1 k/m

Figure 1: Simulink model of a simple oscillator with c/m = 0.1 and k/m = 1 Generating output and plots [t x] = sim(oscillator) plot(t,x(1,:)) Getting state info [sizes xo states] = oscillator sizes = 2 0 0 0 0 0 1 xo = 1 0 states = oscillator/y oscillator/ydot %Runs simulation

Pendulum
+ c + W l sin = u J

Here we consider the motion of a planar pendulum subject to a constant torque u:

Note that this can be rewritten as: W l sin ) . = 1 (u c J

Simple planar pendulum (pend)

thetaddot 1 u 0.5 1/J

1 s

thetadot

1 s

theta Scope

0.4

2*sin(u) gravity

Figure 2: Simulink model of a simple planar pendulum with W l = 2, c = 0.4, J = 2 and u = 1. Finding equilibrium and trim conditions trim(pend) ans = 0.0000 0.5236

Pendulum on cart
Pendulum on cart (pendcart)
Run pendcartpar for parameters ydot f(u) 1 s y 1 s

0 u Mux f(u) thetadot 1 s theta 1 s

Figure 3: Simulink model of pendulum on cart

The motion of the pendulum on a cart can be described by + ml sin 2 = u (M + m) y ml cos ml cos y + ml2 + mlg sin = 0 where M is the mass of the cart, m is the pendulum mass, l is distance from cart to pendulum mass, and g is the gravitational acceleration constant. The variables y and are the cart displacement and the pendulum angle, respectively. and y Solving for yields 2 mg sin cos )/(M + m sin 2 ) y = (u ml sin = (cos u ml sin cos 2 (M + m)g sin cos )/(M l + ml sin 2 ) Simulink model. Top function block (u[1]-m*l*sin(u[3])*u[2]*u[2]-m*g*sin(u[3])*cos(u[3]))/(M+m*sin(u[3])2 ) Bottom function block (cos(u[3])*u[1]-m*l*sin(u[3])*cos(u[3])*u[2]*u[2]-(M+m)*g*sin(u[3])) /(M*l+m*l*sin(u[3])2 ) Specifying parameters. Run this le rst before initial simulations. % %pendcartparam.m % %Set parameters for pendulum on cart example M=1 m=1 l=1 g=1 3

Cannonball
V = g sin V 2 /m V = g cos

and

p = V cos = V sin h

where = SCD /2.

u[1]*cos(u[2])

p 1 s

Dynamics

u[1]*sin(u[2])

1 s h

XY Graph

Figure 4: Simulink model of cannonball

-g*sin(u[2]) -kappa*u[1]^2/m

V 1 s

gamma -g*cos(u[2])/u[1] 1 s

Figure 5: Dynamics subsystem

%cannonballpar.m % %Parameters for cannonball g = 9.81 m= 1 rho=0.3809 Cd=0.4 S=0.1 kappa = rho*S*Cd/2

S-Functions

S-functions are useful when one wants to use equations to describe a dynamical system. Recall the simple pendulum example. Here we use an S-function to describe the pendulum dynamics.
pend2
1 u sfun_pend S-Function

theta Scope

Figure 6: Simulink model of a simple planar pendulum with W l = 2, c = 0.4 and J = 2.

The following Matlab M-le, called sfun pend.m, generates a S-function for the simple pendulum with input u and output . You can use this le as a template for all your S-function M-les. You need only change the lines in the boxes. % sfun pend.m % S-function to describe the dynamics of the simple pendulum % Input is a torque and output is pendulum angle function [sys,x0,str,ts] =sfun pend(t,x,u,flag) % % % % % t is time x is state u is input flag is a calling argument used by Simulink. The value of flag determines what Simulink wants to be executed.

switch flag case 0 % Initialization [sys,x0,str,ts]=mdlInitializeSizes; case 1 % Compute xdot sys=mdlDerivatives(t,x,u); case 2 % Not needed for continuous-time systems

case 3 % Compute output sys = mdlOutputs(t,x,u); case 4 % Not needed for continuous-time systems 6

case 9 end

% Not needed here

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % mdlInitializeSizes %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % function [sys,x0,str,ts]=mdlInitializeSizes % % Create the sizes structure sizes=simsizes; sizes.NumContStates = 2; sizes.NumDiscStates = 0; sizes.NumOutputs = 1; sizes.NumInputs = 1; % Set number of output variables % Set number of input variables % Set number of continuous-time state variables

sizes.DirFeedthrough = 0; sizes.NumSampleTimes = 1; % Need at least one sample time sys = simsizes(sizes); % x0=[0;0]; % Set initial state % str=[]; % str is always an empty matrix ts=[0 0]; % ts must be a matrix of at least one row and two columns % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % mdlDerivatives %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % function sys = mdlDerivatives(t,x,u) % % Compute xdot based on (t,x,u) and set it equal to sys % sys(1) = x(2); sys(2) = (-2*sin(x(1))-0.4*x(2) + u)/2; % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 7

% mdlOutput %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % function sys = mdlOutputs(t,x,u) % % Compute output based on (t,x,u) and set it equal to sys sys = x(1);

Using ODE les in S-functions

sfun_pend3 S-Function

Demux

theta-dot

Run pendpar before initial simulation

Figure 7: Another Simulink model of the simple planar pendulum

The following Matlab M-le generates a S-function for the simple planar pendulum with no . The applied torque u is treated as a parameter. input and two output variables and This example also illustrates how one can use ode les in an S-function.

% sfun_pend3.m % S-function to describe the dynamics of % SIMPLE PLANAR PENDULUM a

%CHANGE

%CHANGE

function [sys,x0,str,ts] = sfun_pend3(t,x,u,flag) % % % % % t is time x is state u is input flag is a calling argument used by Simulink. The value of flag determines what Simulink wants to be executed.

%CHANGE

switch flag case 0 % Initialization [sys,x0,str,ts]=mdlInitializeSizes; case 1 % Compute xdot sys=mdlDerivatives(t,x,u); case 2 case 3 % Not needed for continuous-time systems % Compute output 9

sys = mdlOutputs(t,x,u); case 4 case 9 end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % mdlInitializeSizes %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % function [sys,x0,str,ts]=mdlInitializeSizes % % Create the sizes structure sizes=simsizes; sizes.NumContStates = 2; %Set number of continuous-time state variables %CHANGE sizes.NumDiscStates = 0; sizes.NumOutputs = 2; %Set number of outputs %CHANGE sizes.NumInputs = 0; %Set number of inputs %CHANGE sizes.DirFeedthrough = 0; sizes.NumSampleTimes = 1; sys = simsizes(sizes); % x0=[1;0]; % Not needed for continuous-time systems % Not needed here

%Need at least one sample time

% Set initial state

%CHANGE

str=[]; % str is always an empty matrix ts=[0 0]; % ts must be a matrix of at least one row and two columns % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % mdlDerivatives %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % function sys = mdlDerivatives(t,x,u) % % Compute xdot based on (t,x,u) and set it equal to sys % sys= pendode(t,x); % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 10

%CHANGE

% mdlOutput %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % function sys = mdlOutputs(t,x,u) % % Compute xdot based on (t,x,u) and set it equal to sys sys = [x(1); x(2)]; %CHANGE

The above S-function uses the following ODE le %pendode.m function xdot = pendode(t,x) global W l J c u xdot(1) = x(2); xdot(2) = (-W*l*sin(x(1)) -c*x(2) + u)/J; This ODE le requires the following parameter le %pendpar.m global W l J c u J = 2 c=0.4 W= 2 l=1;

11

Trim and Linearization

Recall that the simple planar pendulum is described by

Hence

For the chosen parameters A= Trim xe=trim(pend3) xe = 0.5236 -0.0000 Linearization A=linmod(pend3,xe) A = 0 1.0000 -0.8660 -0.2000

More u=0; xe=trim(pend3) xe = 1.0e-023 * -0.6617 -0.0009 A=linmod(pend3,xe) A = 0 1.0000 -1.0000 -0.2000

12

Input output stu

theta

1 In1

sfun_pend4 S-Function

1 Out1

Run pend4par before initial simulation

Figure 8: Yet another Simulink model of a simple planar pendulum

Recall that the simple planar pendulum is described by + c + W l sin = u J We will view this as an input output system with input u and output y = . as states. We choose x1 = and x2 = Trim. Suppose that we wish that y = y e = /6 = 0.5236 when the system is in equilibrium. So, we need to compute ue and xe so that y e = 1. From the dierential equation above, we obtain ue = W l sin(y e ) = (2)(1) sin(/6) = 1; also xe = Using the trim command, we obtain [xe ue ye xdote] = trim(pend4,[],[],[pi/6],[],[],[1]) xe = 0.5236 0 ue = ye = xdote = 0 0 1.0000 0.5236 ye 0 = 0.5236 0

13

Linearization. The linearized system is described by x 1 = x2 x 2 = (W l cos y e /J )x1 (c/J )x2 + (1/J )u y = x1 Hence, A= 0 1 W l cos y /J c/J
e

B=

0 1/J

C=

1 0

D = 0.

For the chosen parameters and trim conditions: A= 0 1 0.8660 0.2 , B= 0 1/J , C= 1 0 , D = 0.

MATLAB yields: [A B C D] = linmod(pend4,xe,ue) A = 0 -0.8660 B = 0 0.5000 C = 1.0000 D = 0 lambda =eig(A) lambda = -0.1000 + 0.9252i -0.1000 - 0.9252i Transfer function and poles and zeros. [num den]=ss2tf(A,B,C,D) num = den = 0 1.0000 0 0.2000 0.5000 0.8660 14 0 1.0000 -0.2000

poles=roots(den) poles = -0.1000 + 0.9252i -0.1000 - 0.9252i zeros=roots(num) zeros = Empty matrix: 0-by-1

15