Simulation of control systems in MATLAB/Simulink

Why it matters

Nobody commissions a servo drive or a robot joint straight from hand calculations. Engineers model the plant, close the loop and check step response, margins and actuator effort in simulation first, because it is cheap, fast and safe. MATLAB's Control System Toolbox and Simulink are the standard tools in industry and in university labs (Python's python-control and Scilab/Xcos work the same way), and being able to check a hand result against a simulation is a core interview and lab skill.

Key ideas

Two ways of working.

• MATLAB scripts (Control System Toolbox): build LTI objects and analyse them with commands. Best for linear analysis — poles, step response, Bode, root locus, margins, design.
• Simulink: drag-and-connect block diagrams solved numerically in time. Best when the model has nonlinearities (saturation, dead zone, friction, quantisation), delays, mixed continuous and discrete parts, or several interacting subsystems.

Building models in MATLAB.

• G = tf([1 3],[1 4 5]) creates (s + 3)/(s² + 4s + 5); coefficient vectors are in descending powers of s.
• s = tf('s'); G = 10/(s*(s+2)) — the same thing written algebraically.
• zpk([-3],[-2+1i -2-1i],1) gives zero–pole–gain form; ss(A,B,C,D) gives state space; tf(sys) and ss(sys) convert.
• Interconnection: series(G1,G2) or G1*G2; parallel(G1,G2) or G1+G2; feedback(G,H) for negative feedback, feedback(G,H,+1) for positive.
• Analysis: pole, zero, damp (ωn and ζ of each pole), dcgain, step, impulse, lsim (arbitrary input), stepinfo (rise time, settling time, overshoot), bode, margin (GM, PM and crossovers), nyquist, rlocus, rlocfind.
• Design: pidtune and the PID Tuner app, place and acker for pole placement, ctrb and obsv with rank, c2d(G,T,'zoh') for discretisation.

Building models in Simulink. Typical blocks: Step (source), Sum (set signs "+−" for negative feedback), Gain, Transfer Fcn, State-Space, PID Controller, Integrator, Saturation, Transport Delay, Zero-Order Hold, Scope and To Workspace. Set the Step block's step time and final value, the Sum signs, and the simulation stop time; log signals to compare with hand results.

Solvers. Simulink integrates the model's differential equations numerically. Variable-step solvers (ode45 by default; ode15s for stiff systems with widely separated time constants) adapt the step for accuracy. Fixed-step solvers (ode1 Euler, ode4 Runge–Kutta) are used for code generation and real-time targets. The step size must be small compared with the fastest time constant: explicit Euler on dx/dt = −a·x is stable only for step h < 2/a, and accurate only for h much smaller (about a tenth of the time constant or less).

Good practice. Check every simulation against something you can compute by hand: DC gain, pole locations, final value, overshoot. Model actuator limits (saturation) — a linear model happily asks for 1000 V. Use realistic sensor noise and sample times when testing a digital controller. Keep units consistent throughout (SI).

Formulas

G(s) = num(s)/den(s) — coefficients in descending powers of s, as entered in tf(num, den).

T(s) = G(s)/(1 + G(s)·H(s)) — what feedback(G,H) computes.

ωn = √(den₀), ζ = den₁/(2ωn) for den = s² + den₁·s + den₀ — what damp reports for a complex pair.

DC gain = G(0) = num₀/den₀ — what dcgain reports (stable G).

Explicit Euler: x_(k+1) = x_k + h·f(x_k, u_k); stable for dx/dt = −a·x only if h < 2/a (h in s, a in s⁻¹).

Worked examples

Example 1 (standard — check a step response). Simulate G(s) = 1/(s² + 3s + 2) and predict what MATLAB reports.

1. Script: G = tf(1,[1 3 2]); step(G); stepinfo(G); damp(G); dcgain(G).
2. Poles: s² + 3s + 2 = (s + 1)(s + 2) → −1 and −2, both real, so damp shows ζ = 1 for each pole (overdamped overall, no overshoot).
3. DC gain = 1/2 = 0.5, so the step response settles at 0.5, not 1.
4. Analytic response: y(t) = 0.5 − e^(−t) + 0.5e^(−2t).
5. 10–90 % rise time: y reaches 0.05 at t = 0.380 s and 0.45 at t = 2.970 s → tr ≈ 2.59 s.
6. 2 % settling time (|y − 0.5| ≤ 0.01): ts ≈ 4.60 s; overshoot 0 %. stepinfo should report the same values; if it does not, the model was entered wrongly.
7. Simulink version: Step → Transfer Fcn (numerator [1], denominator [1 3 2]) → Scope, stop time 10 s.

Example 2 (GATE level — closed loop and margins). For G(s) = 10/(s(s + 2)) with unity feedback, predict the output of T = feedback(G,1); damp(T); stepinfo(T); margin(G).

1. T(s) = 10/(s² + 2s + 10): ωn = √10 = 3.16 rad/s, ζ = 2/(2 × 3.16) = 0.316.
2. Overshoot Mp = exp(−π × 0.316/√(1 − 0.1)) = 35.1 %.
3. Peak time tp = π/(ωn√(1 − ζ²)) = π/3.0 = 1.05 s.
4. Settling time (2 % envelope estimate) ≈ 4/(ζωn) = 4/1 = 4 s (stepinfo, which uses the true curve, will give a somewhat smaller value).
5. Margins of G: gain crossover where 10 = ω√(ω² + 4) → ω⁴ + 4ω² − 100 = 0 → ω² = 8.20, ω_gc = 2.86 rad/s; PM = 180° − 90° − tan⁻¹(2.86/2) = 90° − 55.1° = 34.9°. The phase never reaches −180°, so GM = ∞ (margin reports Inf).
6. Steady-state error to a unit ramp: Kv = 10/2 = 5 s⁻¹ → ess = 0.2 — check with lsim(T, t, t).

Common mistakes

• Entering coefficients in ascending order or dropping a zero coefficient: s² + 5 is [1 0 5], not [1 5].
• Using G*H/(1+G*H) by hand instead of feedback, which creates non-minimal models with cancelled poles; use minreal if needed.
• Wrong signs on the Simulink Sum block, giving positive feedback.
• Too short a stop time, so the response has not settled, or too large a fixed step, giving numerical instability that is not physical.
• Trusting a linear simulation that demands impossible actuator signals; add a Saturation block.
• Confusing margin of the open loop (correct) with margin of the closed loop (meaningless).

For GATE ME

GATE does not test software syntax, but simulation-based thinking is exactly what its numerical questions check: ωn, ζ, overshoot, settling time, DC gain, poles and margins of a given loop. Use this topic to practise predicting the result before you simulate, and to verify hand answers for problems from every other topic in this subject.

Quick check

1. What does tf(1,[1 2]) represent, and what is its time constant?
2. What is ωn of 5/(s² + 6s + 5)?
3. Which command gives the closed-loop transfer function for unity negative feedback around G?
4. What is the damping ratio of 1/(s² + 4s + 3)?
5. For explicit Euler on dx/dt = −10x, what step size makes the simulation unstable?

Answers: 1. 1/(s + 2); τ = 0.5 s. 2. √5 = 2.24 rad/s. 3. feedback(G,1). 4. ωn = √3, ζ = 4/(2√3) = 1.15 (overdamped). 5. Any h > 0.2 s.