Problem 459052 · hard · Level 04 Non-Linear Data Structures

How Close Is the Shower to Oscillating?

frequency response · complex numbers · gain margin · phase margin · stability

A thermostatic shower's controller adjusts a mixing valve. Its loop is slow and late: the valve and the pipe warm up with time constants lags (seconds), and the mixed water takes delay seconds to travel to the shower head. Turn the gains up too far and the temperature swings hot and cold by itself. How much room is left is measured by two stability margins, read from the loop's frequency response.

Send a sine wave of angular frequency w (radians per second) round the loop. Every block multiplies it by a complex number: the number's size abs(...) scales the wave and its angle cmath.phase(...) shifts it (negative angles are lags). With s = 1j * w:

controller   C = Kp + Ki / s + Kd * s
plant        K * e^(-s * delay) / ((1 + s * T1) * (1 + s * T2) * ...)   for the T's in lags
loop         L = controller times plant

The size of L is the product of the sizes (the delay has size 1). Its angle is the sum of the angles, which you must add up block by block, in degrees: cmath.phase(C), minus cmath.phase(1 + s * T) for each lag, minus w * delay for the delay (converted from radians). The total can go far below -180 degrees, which cmath.phase of the product cannot show.

Use the grid of n frequencies w[k] = w_lo * (w_hi / w_lo) ** (k / (n - 1)), k = 0 .. n - 1 (evenly spaced on a log scale), and find:

  • the gain crossover: the first k >= 1 with abs(L) at w[k-1] above 1 and at w[k] at most 1. Its phase margin is 180 + angle(L at w[k]) degrees: how much more lag the loop could take before a returning wave arrives exactly upside down and feeds itself;
  • the phase crossover: the first k >= 1 with the angle at w[k-1] above -180 and at w[k] at most -180. Its gain margin is -20 * log10(abs(L at w[k])) dB: how much the gain could grow before the loop oscillates.

Write margins(Kp, Ki, Kd, K, lags, delay, w_lo, w_hi, n) that returns the tuple (w_gain, phase_margin, w_phase, gain_margin), with None for both values of a crossover that does not happen on the grid. Press Run with plot(ws, sizes, logx=True) and the angles to draw the loop's Bode plot.

Examples

Input:  Kp = 0.8, Ki = 0.2, Kd = 0, K = 1.5, lags = [4, 1], delay = 0.8,
        w_lo = 0.01, w_hi = 10, n = 400
Output: (0.29252663, 60.28596307, 1.0, 13.46787486)
Explanation: the loop gain falls through 1 near 0.29 rad/s, where the lag is 119.7 degrees:
60.3 degrees of margin. The lag reaches 180 degrees near 1 rad/s, where the loop gain is
0.21, so the gain could grow 13.5 dB (4.7 times) before the shower oscillates.

Input:  Kp = 1, Ki = 0, Kd = 0, K = 3, lags = [], delay = 0.5, w_lo = 0.1, w_hi = 100, n = 50
Output: (None, None, 6.86648845, -9.54242509)
Explanation: a pure delay with a loop gain of 3 never drops below 1, and at the phase
crossover the gain margin is negative: this loop oscillates.

Constraints

  • Kp >= 0, K > 0, w_lo < w_hi, and the controller is not all zero (Kp + Ki + Kd > 0)
  • answers are compared with a tolerance of 1e-6; do not round

Goals

  • Compute a loop's frequency response with Python's complex numbers
  • Add the phases of the blocks in a loop instead of taking the phase of their product
  • Find the gain and phase crossovers on a frequency grid and read off both margins
Starting Python…