Merge pull request #1 from edrivedesign/copilot/create-control-system-model

edrivedesign · web-flow · commit 170dfed8e459 · 2026-07-01T20:26:22.000+01:00
Add steam engine governor simulation example with interactive GUI
diff --git a/examples/steam_engine_governor.py b/examples/steam_engine_governor.py
@@ -0,0 +1,310 @@
+#!/usr/bin/env python3
+"""Steam engine governor model with a simple input/output GUI."""
+
+import os
+
+if 'PYCONTROL_TEST_EXAMPLES' in os.environ:
+    import matplotlib
+    matplotlib.use('Agg')
+
+import control as ct
+import matplotlib.pyplot as plt
+import numpy as np
+
+try:
+    import tkinter as tk
+    from tkinter import ttk
+    from matplotlib.backends.backend_tkagg import FigureCanvasTkAgg
+except ImportError:
+    tk = None
+    ttk = None
+    FigureCanvasTkAgg = None
+
+
+DEFAULTS = {
+    'kp': 1.0,
+    'ti': 2.0,
+    'tg': 0.1,
+    'ts': 0.2,
+    'kt': 2.0,
+    'j': 4.5,
+    'b': 1.2,
+    'ref_step': 0.08,
+    'load_step': 0.10,
+    'load_time': 5.0,
+    'sim_time': 25.0,
+}
+SETTLING_TOLERANCE = 0.02
+MIN_SETTLING_TOLERANCE = 1e-3
+NUM_DOMINANT_POLES = 2
+
+
+def create_governor_system(params):
+    """Create a linearized steam engine governor model."""
+    kp = params['kp']
+    ti = params['ti']
+    tg = params['tg']
+    ts = params['ts']
+    kt = params['kt']
+    inertia = params['j']
+    damping = params['b']
+    for name, value in (
+        ('ti', ti), ('tg', tg), ('ts', ts), ('kt', kt), ('j', inertia)):
+        if value <= 0:
+            raise ValueError(f'{name} must be positive')
+    if damping < 0:
+        raise ValueError('b must be nonnegative')
+    ki = kp / ti
+
+    a_matrix = np.array([
+        [0, 0, 0, -1],
+        [ki / tg, -1 / tg, 0, -kp / tg],
+        [0, kt / ts, -1 / ts, 0],
+        [0, 0, 1 / inertia, -damping / inertia],
+    ])
+    b_matrix = np.array([
+        [1, 0],
+        [kp / tg, 0],
+        [0, 0],
+        [0, -1 / inertia],
+    ])
+    c_matrix = np.array([
+        [0, 0, 0, 1],
+        [0, 1, 0, 0],
+        [0, 0, 1, 0],
+        [0, 0, 0, -1],
+    ])
+    d_matrix = np.array([
+        [0, 0],
+        [0, 0],
+        [0, 0],
+        [1, 0],
+    ])
+
+    return ct.ss(
+        a_matrix, b_matrix, c_matrix, d_matrix,
+        states=['integral_error', 'valve', 'steam_torque', 'speed'],
+        inputs=['reference', 'load_torque'],
+        outputs=['speed', 'valve', 'steam_torque', 'error'],
+        name='steam_governor',
+    )
+
+
+def simulate_governor(params, samples=600):
+    """Simulate the governor response for the requested settings."""
+    system = create_governor_system(params)
+    time = np.linspace(0, params['sim_time'], samples)
+    reference = np.full_like(time, params['ref_step'])
+    load = np.where(time >= params['load_time'], params['load_step'], 0.0)
+    response = ct.input_output_response(system, time, [reference, load])
+    return system, time, np.asarray(response.outputs), reference, load
+
+
+def settling_time(time, signal, target, tolerance):
+    """Return the first time after which the response stays in band."""
+    error = np.abs(signal - target)
+    indices = np.nonzero(error > tolerance)[0]
+    if indices.size == 0:
+        return time[0]
+    last_outside = indices[-1]
+    if last_outside >= len(time) - 1:
+        return np.nan
+    return time[last_outside + 1]
+
+
+def summarize_response(system, time, outputs, target):
+    """Create summary metrics for the speed response."""
+    speed = outputs[0]
+    final_speed = speed[-1]
+    peak_speed = np.max(speed)
+    minimum_speed = np.min(speed)
+    tolerance = max(
+        SETTLING_TOLERANCE * max(abs(target), abs(final_speed)),
+        MIN_SETTLING_TOLERANCE)
+    if abs(target) < 1e-9:
+        overshoot = 0.0
+    elif target > 0:
+        overshoot = max(0.0, (peak_speed - target) / abs(target) * 100)
+    else:
+        overshoot = max(0.0, (target - minimum_speed) / abs(target) * 100)
+    steady_state_error = target - final_speed
+    poles = list(ct.poles(system))
+    # Use the least stable (rightmost) poles as the dominant dynamics.
+    dominant_poles = sorted(
+        poles, key=lambda pole: pole.real)[-min(NUM_DOMINANT_POLES, len(poles)):]
+    damping_ratios = []
+    for pole in dominant_poles:
+        magnitude = abs(pole)
+        if magnitude > 0:
+            damping_ratios.append(max(0.0, min(1.0, -pole.real / magnitude)))
+
+    return {
+        'settling_time': settling_time(time, speed, target, tolerance),
+        'overshoot': overshoot,
+        'steady_state_error': steady_state_error,
+        'dominant_poles': dominant_poles,
+        'damping_ratios': damping_ratios,
+    }
+
+
+def summary_text(metrics):
+    """Format the response metrics for display."""
+    poles = ', '.join(
+        f'{pole.real:.2f}{pole.imag:+.2f}j' for pole in metrics['dominant_poles'])
+    if metrics['damping_ratios']:
+        damping = ', '.join(f'{ratio:.2f}' for ratio in metrics['damping_ratios'])
+    else:
+        damping = 'n/a'
+    settling = (
+        f"{metrics['settling_time']:.2f} s"
+        if np.isfinite(metrics['settling_time']) else 'not settled'
+    )
+    return '\n'.join([
+        f'Settling time ({SETTLING_TOLERANCE:.0%} band): {settling}',
+        f'Overshoot: {metrics["overshoot"]:.2f}%',
+        f'Steady-state error: {metrics["steady_state_error"]:.4f}',
+        f'Dominant poles: {poles}',
+        f'Damping ratios: {damping}',
+    ])
+
+
+def plot_results(figure, time, outputs, reference, load, metrics):
+    """Draw the response plots for the governor simulation."""
+    speed, valve, torque, error = outputs
+    figure.clf()
+    axes = figure.subplots(3, 1, sharex=True)
+
+    axes[0].plot(time, speed, label='Speed deviation')
+    axes[0].plot(time, reference, '--', label='Reference change')
+    axes[0].set_ylabel('Speed')
+    axes[0].set_title('Steam engine governor response')
+    axes[0].grid(True, linestyle='dotted')
+    axes[0].legend(loc='best')
+
+    axes[1].plot(time, valve, label='Valve position')
+    axes[1].plot(time, torque, label='Steam torque')
+    axes[1].set_ylabel('Actuation')
+    axes[1].grid(True, linestyle='dotted')
+    axes[1].legend(loc='best')
+
+    axes[2].plot(time, load, label='Load torque disturbance')
+    axes[2].plot(time, error, '--', label='Speed error')
+    axes[2].set_xlabel('Time [s]')
+    axes[2].set_ylabel('Disturbance / error')
+    axes[2].grid(True, linestyle='dotted')
+    axes[2].legend(loc='best')
+
+    figure.text(
+        0.99, 0.98, summary_text(metrics), va='top', ha='right',
+        bbox={'boxstyle': 'round', 'facecolor': 'white', 'alpha': 0.9},
+    )
+    figure.tight_layout(rect=[0, 0, 0.86, 1])
+
+
+class GovernorApp:
+    """GUI for adjusting governor inputs and viewing outputs."""
+
+    field_specs = (
+        ('Governor proportional gain (Kp)', 'kp'),
+        ('Integral time (Ti) [s]', 'ti'),
+        ('Actuator time constant (Tg) [s]', 'tg'),
+        ('Steam lag (Ts) [s]', 'ts'),
+        ('Steam torque gain (Kt)', 'kt'),
+        ('Equivalent inertia (J)', 'j'),
+        ('Damping coefficient (B)', 'b'),
+        ('Reference step', 'ref_step'),
+        ('Load torque step', 'load_step'),
+        ('Load step time [s]', 'load_time'),
+        ('Simulation time [s]', 'sim_time'),
+    )
+
+    def __init__(self, root):
+        self.root = root
+        self.root.title('Steam Engine Governor')
+        self.variables = {
+            key: tk.StringVar(value=str(value)) for key, value in DEFAULTS.items()
+        }
+        self.status = tk.StringVar(
+            value='Adjust the model inputs and click Simulate.')
+
+        main = ttk.Frame(root, padding=10)
+        main.pack(fill='both', expand=True)
+        controls = ttk.LabelFrame(main, text='Model inputs', padding=10)
+        controls.pack(side='left', fill='y')
+        plots = ttk.Frame(main)
+        plots.pack(side='right', fill='both', expand=True)
+
+        for row, (label, key) in enumerate(self.field_specs):
+            ttk.Label(controls, text=label).grid(
+                row=row, column=0, sticky='w', padx=(0, 8), pady=3)
+            ttk.Entry(controls, textvariable=self.variables[key], width=12).grid(
+                row=row, column=1, sticky='ew', pady=3)
+
+        ttk.Button(controls, text='Simulate', command=self.update_plots).grid(
+            row=len(self.field_specs), column=0, columnspan=2, sticky='ew',
+            pady=(10, 4))
+        ttk.Label(
+            controls, textvariable=self.status, justify='left', wraplength=260
+        ).grid(row=len(self.field_specs) + 1, column=0, columnspan=2, sticky='w')
+
+        self.figure = plt.Figure(figsize=(10, 7))
+        self.canvas = FigureCanvasTkAgg(self.figure, master=plots)
+        self.canvas.get_tk_widget().pack(fill='both', expand=True)
+
+        self.update_plots()
+
+    def read_params(self):
+        """Read and validate values from the input fields."""
+        params = {key: float(variable.get()) for key, variable in self.variables.items()}
+        for key in ('ti', 'tg', 'ts', 'kt', 'j', 'sim_time'):
+            if params[key] <= 0:
+                raise ValueError(f'{key} must be positive')
+        if params['b'] < 0:
+            raise ValueError('b must be nonnegative')
+        if params['load_time'] < 0:
+            raise ValueError('load_time must be nonnegative')
+        if params['load_time'] > params['sim_time']:
+            raise ValueError('load_time must not exceed sim_time')
+        return params
+
+    def update_plots(self):
+        """Recompute the simulation and refresh the output plots."""
+        try:
+            system, time, outputs, reference, load = simulate_governor(
+                self.read_params())
+            metrics = summarize_response(system, time, outputs, reference[-1])
+        except Exception as exc:
+            self.status.set(f'Unable to simulate: {exc}')
+            return
+
+        plot_results(self.figure, time, outputs, reference, load, metrics)
+        self.canvas.draw()
+        self.status.set(summary_text(metrics))
+
+
+def run_headless_demo():
+    """Exercise the model and plotting code without starting the GUI."""
+    figure = plt.figure(figsize=(10, 7))
+    system, time, outputs, reference, load = simulate_governor(DEFAULTS)
+    metrics = summarize_response(system, time, outputs, reference[-1])
+    plot_results(figure, time, outputs, reference, load, metrics)
+    figure.canvas.draw()
+    plt.close(figure)
+
+
+def main():
+    """Start the GUI unless the example test harness is active."""
+    if 'PYCONTROL_TEST_EXAMPLES' in os.environ:
+        run_headless_demo()
+        return
+
+    if tk is None or FigureCanvasTkAgg is None:
+        raise SystemExit('Tkinter support is required to run this GUI example.')
+
+    root = tk.Tk()
+    GovernorApp(root)
+    root.mainloop()
+
+
+if __name__ == '__main__':
+    main()
