import numpy as np import matplotlib.pyplot as plt from scipy.signal import TransferFunction, lsim # ============================================================ # Circuit Parameters # ============================================================ R = 1.5 # Resistance [Ω] L = 0.0025 # Inductance [H] (2.5 mH) V_SUPPLY = 12.0 # Supply rail [V] tau = L / R # RL time constant ≈ 1.667 ms # ============================================================ # PWM Parameters # ============================================================ F_PWM = 16e3 # PWM frequency [Hz] T_PWM = 1.0 / F_PWM # PWM period [s] (62.5 µs) # ============================================================ # Simulation Window # ============================================================ T_END = 1e-3 # 1 ms → 16 full PWM cycles DT = 1e-7 # 100 ns resolution (625 pts / PWM cycle) t = np.arange(0, T_END + DT / 2, DT) # ============================================================ # Duty-Cycle Command D(t) # ============================================================ # Ramp from 20 % → 80 % over the window so every PWM cycle # has a visibly different pulse width. def duty_command(t_val): """Continuous duty-cycle setpoint (from a controller).""" return np.clip(0.2 + 0.6 * (np.asarray(t_val) / T_END), 0.0, 1.0) D_continuous = duty_command(t) # ============================================================ # MODEL 1 — Abstracted (Average-Voltage) Approximation # ============================================================ # Treats the coil voltage as the smooth signal D(t)·V. # Transfer function: I(s)/D(s) = V / (Ls + R) G = TransferFunction([V_SUPPLY], [L, R]) _, i_avg, _ = lsim(G, U=D_continuous, T=t) # ============================================================ # MODEL 2 — True PWM Waveform (exact analytical solution) # ============================================================ # Between every switching edge the RL circuit obeys: # # di/dt = (V_seg − R·i) / L (V_seg = V_SUPPLY or 0) # # Closed-form from initial condition i₀ at time t₀: # # i(t) = V_seg/R + (i₀ − V_seg/R) · exp(−R·(t − t₀) / L) # # We propagate i analytically from edge to edge — zero # numerical-integration error. The only error source is # IEEE-754 floating-point arithmetic (~1e-15 relative). # --- Step 1: build segment table and propagate boundary currents --- seg_t_start = [] # start time of each constant-V segment seg_V = [] # voltage applied during segment seg_i0 = [] # current at segment start i_boundary = 0.0 # coil starts de-energised cycle = 0 while cycle * T_PWM < T_END: t_cycle = cycle * T_PWM D = float(duty_command(t_cycle)) # ---- ON phase (V_SUPPLY applied) ---- t_on_start = t_cycle t_on_end = min(t_cycle + D * T_PWM, T_END) if t_on_end > t_on_start: seg_t_start.append(t_on_start) seg_V.append(V_SUPPLY) seg_i0.append(i_boundary) dt_on = t_on_end - t_on_start i_boundary = (V_SUPPLY / R) + (i_boundary - V_SUPPLY / R) * np.exp(-R * dt_on / L) # ---- OFF phase (0 V applied, free-wheeling through diode) ---- t_off_start = t_on_end t_off_end = min((cycle + 1) * T_PWM, T_END) if t_off_end > t_off_start: seg_t_start.append(t_off_start) seg_V.append(0.0) seg_i0.append(i_boundary) dt_off = t_off_end - t_off_start i_boundary = i_boundary * np.exp(-R * dt_off / L) cycle += 1 seg_t_start = np.array(seg_t_start) seg_V = np.array(seg_V) seg_i0 = np.array(seg_i0) # --- Step 2: evaluate on the dense time array (vectorised) --- idx = np.searchsorted(seg_t_start, t, side='right') - 1 idx = np.clip(idx, 0, len(seg_t_start) - 1) dt_in_seg = t - seg_t_start[idx] V_at_t = seg_V[idx] i0_at_t = seg_i0[idx] i_pwm = (V_at_t / R) + (i0_at_t - V_at_t / R) * np.exp(-R * dt_in_seg / L) v_pwm = V_at_t # switching waveform for plotting v_avg = D_continuous * V_SUPPLY # average-model voltage # ============================================================ # Console Output — sanity-check steady-state values # ============================================================ print(f"RL time constant τ = L/R = {tau*1e3:.3f} ms") print(f"PWM period T = 1/f = {T_PWM*1e6:.1f} µs") print(f"Sim resolution Δt = {DT*1e9:.0f} ns ({int(T_PWM/DT)} pts/cycle)") print() print("Expected steady-state currents i_ss = (V/R)·D :") for d in [0.2, 0.4, 0.6, 0.8]: print(f" D = {d:.1f} → i_ss = {V_SUPPLY / R * d:.3f} A") # ============================================================ # Plotting — 4-panel comparison # ============================================================ t_us = t * 1e6 # time axis in µs fig, axes = plt.subplots(4, 1, figsize=(14, 10), sharex=True) fig.suptitle("PWM RL-Circuit Model Comparison (16 kHz, 1 ms window)", fontsize=13, fontweight='bold') # --- 1. Duty-cycle command --- ax = axes[0] ax.plot(t_us, D_continuous * 100, color='tab:purple', linewidth=1.2) ax.set_ylabel("Duty Cycle [%]") ax.set_ylim(0, 100) ax.grid(True, alpha=0.3) # --- 2. Voltage waveforms --- ax = axes[1] ax.plot(t_us, v_pwm, color='tab:orange', linewidth=0.6, label="True PWM voltage") ax.plot(t_us, v_avg, color='tab:blue', linewidth=1.4, label="Average voltage D·V", linestyle='--') ax.set_ylabel("Voltage [V]") ax.set_ylim(-0.5, V_SUPPLY + 1) ax.legend(loc='upper left', fontsize=9) ax.grid(True, alpha=0.3) # --- 3. Current comparison --- ax = axes[2] ax.plot(t_us, i_pwm, color='tab:red', linewidth=0.7, label="True PWM current (exact)") ax.plot(t_us, i_avg, color='tab:blue', linewidth=1.4, label="Averaged-model current", linestyle='--') ax.set_ylabel("Current [A]") ax.legend(loc='upper left', fontsize=9) ax.grid(True, alpha=0.3) # --- 4. Difference / ripple --- ax = axes[3] ax.plot(t_us, (i_pwm - i_avg) * 1e3, color='tab:green', linewidth=0.7) ax.set_ylabel("Δi (PWM − avg) [mA]") ax.set_xlabel("Time [µs]") ax.grid(True, alpha=0.3) plt.tight_layout() plt.show()