MATLAB Training 2: Systems Modeling

Originally copied from a Markdown renderer with LaTeX support

Systems Modeling

Basic steps to systems modeling:

Find the system's coordinates (positions, voltages, pressures, etc.)
Describe the force laws governing the derivatives of the coordinates
Join those force laws into equations of motion (EOM) using the equation of state
- $\sum F = ma$ (mechanical-translational)
- Mass conservation (fluid, pneumatic)
- KVL (electrical)
Numerically integrate the equations of motion to find the trajectory
! SANITY CHECK RESULTS !

Example 1: Spring-mass-damper

Linear spring: force is proportional to displacement
Linear damper: force is proportional to displacement velocity
NOTE: Currently position is displacement, but if there are two coupled
bodies, it doesn't have to be

Translating spring-mass-damper

Equation of motion

From Newton's first law:

$$\sum F = m\ddot{x} = -b\dot{x} - kx + f$$

dt = 0.01;
t_f = 10;
t = 0:dt:t_f;

m = 2; % [kg] mass
k = 20; % [N/m] spring constant
b = 5; % [N/(m/s)] damping constant

f_step = 10 * heaviside(t - 5); % [N] applied force samples
x_0 = 2;
x_dot_0 = -10;

EOM Euler

x_step = zeros(size(t));
x_dot_step = zeros(size(t));

x_step(1) = x_0;
x_dot_step(1) = x_dot_0;

for i = 2:length(t)
    x_prev = x_step(i-1);
    x_dot_prev = x_dot_step(i-1);
    f_prev = f_step(i-1);

    a_t = (-b*x_dot_prev - k*x_prev + f_prev)/m;

    x_dot_step(i) = x_dot_prev + a_t * dt;
    x_step(i) = x_prev + x_dot_prev * dt;
end

Results and sanity checks

Common pitfalls: Forgetting to multiply/divide, putting wrong sign on forces

Stable or unstable?
- Stable
Steady-state condition: Constant, what is the static gain?
- Use physical intuition: Not moving, not accelerating
- Balance forces, calculate using spring constant
Undamped, overdamped, critically damped?
- Calculate the damping ratio from the ODE's characteristic equation
- Only useful for linear second-order systems

Simulink block diagram

Block diagram basics

Generally:

Centered around the integrator
You have a formula for the highest derivative, so this allows you to get
all the lower-order derivatives
Sum, product, gain work as you expect

Simulink block diagrams:

Take variable names from the MATLAB workspace
Use "To Workspace" blocks to put data back to MATLAB

State selection

The equations of motion will be in terms of the system’s states - the variables that store energy. In this case, we have two: position and velocity. We then express (in signals) their derivatives. Velocity’s is acceleration, a non-trivial function of the states and some parameters. Position’s is trivially velocity.

Building the diagram

Begin with the series of integrators, representing the system states. Express the first integrator’s input (acceleration) in terms of the forces, which in turn rely on the outputs of the integrators. The second integrator’s input is just the first integrator’s output, becuase position and velocity are trivially related.

Translating spring-mass-damper block diagram

It looks like you’ve drawn circles, and it’s not necessarily clear how Simulink figures out what to do. In a “visual programming langauge” (LABView) this would be true, but Simulink solves this with blocks that have initial conditions. The one we care about is the integrator, which breaks the algebraic loop by initializing its output to the IC. Then, all of the signals and the integrators' inputs are calculated, the integrators' new state is calculated, and the cycle begins anew.

Note on integrator ICs

Integrators need initial conditions (more on that later). There are two ways of providing them, internally and externally. Internally, you double-click to set. Externally, a signal provides this value, and this signal is read every time the integrator is initialized (this can happen at runtime if you configure integrator reset conditions).

I don’t like internal ICs because it’s easy to forget to set them, and when you look at a block diagram, it is not immediately apparent that you did.

Providing the signal externally (through a Constant block that accesses a workspace variable) makes it clear whether or not the modeller remembered to correctly initialize the integrator.

Plumbing

Write the variables into your MATLAB workspace
Inputs should be timeseries to avoid any ambiguity
Set up the solver with t_f end time, fixed-step ode1 solver, dt step
Use sim() with the block diagram's name
Output is a structure
- A few defaults tout, yout, etc.
- "To workspace" blocks saved as their variable names

force = timeseries(f_step, t);
smd_out = sim("smd.slx");
figure("name", "Example 1 Simulink output");
plot(smd_out.x);

Sanity checks

Same as Euler?
Yes
What if we decrese the damping?
What if we increase the spring stiffness?
What if we double the applied force?

Summary

Revisiting the steps:

Coordinate: $x$
Force laws come from the spring $kx$, damper $b\dot{x}$
Equation of state $\sum F = ma$ gives $m\ddot{x} + b\dot{x} + kx = f$
Euler for-loop, Simulink B.D. to integrate
System has correct stability, steady-state, damping

Example 2: Rotary spring-mass-damper

Spinning this time
- What a twist
Different coordinate systems: displacment and angle
- $x = r\theta$
- $v = r\omega$
- (under a small-angle approximation)
Equation of state: $\sum M = J\ddot{\theta}$

Equation of motion

Spring displacement $L\theta$
Damper speed $(L+l)\theta$
Join using $\sum M = J\ddot{\theta}$

$$ J\ddot{\theta} + b(L+l)^2\dot{\theta} + kL^2\theta = -f(L+l)$$

Block diagram

Summary

Coordinate: $\theta$ (But a transformation (helped by small-angle assumption) to use
translational elements)
Force laws come from the spring $kx$, damper $b\dot{x}$
Equation of state $\sum F = ma$ gives $m\ddot{x} + b\dot{x} + kx = f$
Euler for-loop, Simulink B.D. to integrate
System has correct stability, steady-state, damping

Example 3: Coupled spring-mass-damper

Equations of motion

$k_2$ displacement is $L\theta - x$
$b_2$ displacement is $(L+l)\dot{\theta} - \dot{x}$
Find the correct sign convention by thinking about what happens when one
coordinate is zero

Translating mass:

$$m\ddot{x} = -k_1 x - b_1 \dot{x} - k_2(L\theta - x) - b_2((L+l)\dot\theta - \dot{x})$$

Rotating inertia:

$$J\ddot{\theta} = -k_2 (L\theta - x) - b_2 ((L+l)\dot\theta - \dot{x}) - (L+l)f$$

Block diagram

Summary

Coordinates: $x$ and $\theta$
Force laws come from the springs and dampers
Make sure to use extension instead of position for springs and dampers
Equations of state $\sum F = m\ddot{x}$, $\sum M = J\ddot{\theta}$
Simulink B.D. joining the two models
System has correct stability, steady-state, damping

Summary Summary

You should now know how to:

Come up with equations of motion
How to put EOM into MATLAB
- Especially with coupled states, Simulink is much easier