Engineering Systems

Engineering Systems

Overview

This skill covers computational engineering: control systems design and analysis, signal processing, reliability engineering, mathematical optimization, and sensor data processing. Use it for engineering research problems involving dynamic systems, sensor streams, or system optimization.

When to Use This Skill

• Analyzing or designing control systems (PID, feedback loops)
• Processing sensor time series (filtering, calibration, anomaly detection)
• Running engineering optimization (LP, MIP, scheduling)
• Assessing system reliability or failure analysis
• Working with FEA simulation outputs

1. Control Systems

Transfer Functions and System Analysis

import control
import numpy as np
import matplotlib.pyplot as plt

# Define transfer function G(s) = 1 / (s² + 2s + 1)
G = control.tf([1], [1, 2, 1])
print(G)

# Step response
t, y = control.step_response(G)
rise_time = t[np.argmax(y >= 0.1 * y[-1])]  # 10% rise time (approx)
settling_idx = np.where(np.abs(y - y[-1]) > 0.02 * y[-1])[0]
settling_time = t[settling_idx[-1]] if len(settling_idx) > 0 else 0

print(f"Steady-state value: {y[-1]:.3f}")
print(f"Rise time (10%): {rise_time:.3f} s")
print(f"Settling time (2%): {settling_time:.3f} s")
print(f"Overshoot: {(y.max() - y[-1]) / y[-1] * 100:.1f}%")

# Bode plot: gain margin and phase margin
gm, pm, wcg, wcp = control.margin(G)
print(f"Gain margin: {20*np.log10(gm):.1f} dB  (should be > 6 dB for stability)")
print(f"Phase margin: {pm:.1f}°  (should be > 30° for robustness)")

# Poles and zeros (stability: all poles must have Re < 0 for stability)
poles = control.poles(G)
zeros = control.zeros(G)
stable = all(p.real < 0 for p in poles)
print(f"Poles: {poles}  |  Stable: {stable}")

PID Controller Design

# PID controller: C(s) = Kp + Ki/s + Kd*s
# Closed-loop system: T(s) = C(s)*G(s) / (1 + C(s)*G(s))

def tune_pid_ziegler_nichols(Ku: float, Tu: float) -> dict:
    """Ziegler-Nichols tuning from ultimate gain Ku and period Tu.
    Ku: ultimate gain (gain at stability boundary)
    Tu: ultimate period (oscillation period at Ku)
    """
    return {
        "P":   {"Kp": 0.5 * Ku,      "Ki": 0,                 "Kd": 0},
        "PI":  {"Kp": 0.45 * Ku,     "Ki": 0.54 * Ku / Tu,    "Kd": 0},
        "PID": {"Kp": 0.6 * Ku,      "Ki": 1.2 * Ku / Tu,     "Kd": 3 * Ku * Tu / 40},
    }

# Build and simulate closed-loop PID
def make_pid_controller(Kp, Ki, Kd):
    """Return PID transfer function."""
    # C(s) = Kd*s² + Kp*s + Ki) / s
    return control.tf([Kd, Kp, Ki], [1, 0])

G_plant = control.tf([1], [1, 3, 2])  # Example plant
params = tune_ziegler_nichols = {"Kp": 3.0, "Ki": 1.5, "Kd": 0.5}  # Example values
C = make_pid_controller(**params)
T_cl = control.feedback(C * G_plant, 1)  # Closed-loop
t_cl, y_cl = control.step_response(T_cl)

State-Space Representation

import numpy as np

# dx/dt = Ax + Bu
# y = Cx + Du
A = np.array([[0, 1], [-2, -3]])  # system matrix
B = np.array([[0], [1]])           # input matrix
C = np.array([[1, 0]])             # output matrix
D = np.array([[0]])                # feedthrough

sys = control.ss(A, B, C, D)
print(f"Eigenvalues (poles): {np.linalg.eigvals(A)}")
print(f"Controllable: {np.linalg.matrix_rank(control.ctrb(A, B)) == A.shape[0]}")
print(f"Observable: {np.linalg.matrix_rank(control.obsv(A, C)) == A.shape[0]}")

2. Signal Processing

from scipy import signal
import numpy as np

fs = 1000  # Hz
t = np.arange(0, 5, 1/fs)
# Signal: 10 Hz + 120 Hz (noise) + random noise
raw = (np.sin(2*np.pi*10*t) + 0.3*np.sin(2*np.pi*120*t)
       + 0.5*np.random.randn(len(t)))

# ── Filter Design ─────────────────────────────────────────────
# Butterworth lowpass (maximally flat passband)
sos_lp = signal.butter(N=4, Wn=50, btype="low", fs=fs, output="sos")
filtered_lp = signal.sosfiltfilt(sos_lp, raw)  # zero-phase

# Bandpass filter (8–12 Hz, alpha band example)
sos_bp = signal.butter(N=4, Wn=[8, 12], btype="band", fs=fs, output="sos")
filtered_bp = signal.sosfiltfilt(sos_bp, raw)

# ── FFT Analysis ─────────────────────────────────────────────
N = len(raw)
freqs = np.fft.rfftfreq(N, 1/fs)
fft_mag = np.abs(np.fft.rfft(raw * np.hanning(N)))  # Hanning window

# Find dominant frequencies
peak_idx = np.argsort(fft_mag)[-5:][::-1]
print("Top 5 frequencies:")
for i in peak_idx:
    print(f"  {freqs[i]:.1f} Hz: amplitude {fft_mag[i]:.2f}")

# ── Power Spectral Density (Welch) ────────────────────────────
freq_psd, psd = signal.welch(raw, fs=fs, nperseg=256, noverlap=128)
# More reliable than single FFT; averages over segments

# ── Spectrogram (time-frequency) ─────────────────────────────
freq_s, time_s, Sxx = signal.spectrogram(raw, fs=fs, nperseg=128, noverlap=64)

Filter design rules:

• Always use sosfiltfilt (not lfilter) for zero-phase filtering (no temporal distortion)
• Apply sosfilt when real-time processing is needed (causal, but introduces phase lag)
• Butterworth: flat passband, gentle rolloff — good default
• Chebyshev: steeper rolloff, ripple in passband/stopband — when sharp cutoff matters

3. Reliability Engineering

import math
import numpy as np
from scipy.stats import weibull_min
import matplotlib.pyplot as plt

# Weibull distribution for failure time analysis
# Parameters: shape β (failure mode) and scale η (characteristic life)
# β < 1: decreasing failure rate (infant mortality, defects)
# β = 1: constant failure rate (random failures — exponential distribution)
# β > 1: increasing failure rate (wear-out)

failure_times = np.array([120, 240, 350, 180, 420, 95, 310, 280, 150, 390])  # hours

# Fit Weibull distribution
beta, _, eta = weibull_min.fit(failure_times, floc=0)
print(f"Shape (β): {beta:.2f}  ({'wear-out' if beta>1 else 'early failure' if beta<1 else 'random'})")
print(f"Scale (η): {eta:.1f} hours (characteristic life)")
print(f"MTBF = η·Γ(1+1/β) = {eta * math.gamma(1 + 1/beta):.1f} hours")

# Reliability function R(t) = exp(-(t/η)^β) = P(failure > t)
def reliability(t, beta, eta):
    return np.exp(-(t / eta) ** beta)

t_range = np.linspace(0, 600, 300)
R = reliability(t_range, beta, eta)

# B10 life: time at which 10% of units have failed
from scipy.optimize import brentq
b10 = brentq(lambda t: reliability(t, beta, eta) - 0.90, 0.1, 10000)
print(f"B10 life (10% failure): {b10:.1f} hours")

Fault Tree Analysis (Logic)

# AND gate: system fails if ALL subsystems fail
# OR gate: system fails if ANY subsystem fails

def system_reliability_and(*component_reliabilities):
    """AND gate: all must work for system to work."""
    result = 1.0
    for r in component_reliabilities:
        result *= r
    return result

def system_reliability_or(*component_reliabilities):
    """OR gate: at least one must fail for system to fail (redundant system)."""
    failure_prob = 1.0
    for r in component_reliabilities:
        failure_prob *= (1 - r)
    return 1 - failure_prob  # system reliability

# Example: 2-out-of-3 redundant system
R_comp = 0.95
R_redundant = system_reliability_or(R_comp, R_comp, R_comp)
print(f"2-of-3 redundant system reliability: {R_redundant:.4f}")

4. Engineering Optimization (LP / MIP)

import pulp

# Production planning LP
# Maximize: 25x + 30y (profit)
# Subject to: 2x + y ≤ 100 (labor hours), x + 2y ≤ 80 (materials), x,y ≥ 0

prob = pulp.LpProblem("Production Planning", pulp.LpMaximize)
x = pulp.LpVariable("Product_A", lowBound=0)
y = pulp.LpVariable("Product_B", lowBound=0)

prob += 25*x + 30*y              # objective
prob += 2*x + y <= 100           # labor constraint
prob += x + 2*y <= 80            # materials constraint

prob.solve(pulp.PULP_CBC_CMD(msg=0))
print(f"Status: {pulp.LpStatus[prob.status]}")
print(f"Product A: {x.value():.1f}, Product B: {y.value():.1f}")
print(f"Profit: ${prob.objective.value():.2f}")

# For MIP (integer variables):
# x_int = pulp.LpVariable("x", lowBound=0, cat="Integer")
# y_bin = pulp.LpVariable("y", cat="Binary")  # 0/1 decision

# Google OR-Tools for more complex scheduling/routing:
# from ortools.sat.python import cp_model

5. Sensor Data Processing

import numpy as np
import pandas as pd
from scipy import stats, signal

# ── Calibration ───────────────────────────────────────────────
def linear_calibrate(raw_values: np.ndarray, cal_points: list) -> np.ndarray:
    """Two-point calibration: cal_points = [(raw1, true1), (raw2, true2)]."""
    raw1, true1 = cal_points[0]
    raw2, true2 = cal_points[1]
    slope = (true2 - true1) / (raw2 - raw1)
    intercept = true1 - slope * raw1
    return slope * raw_values + intercept

# ── Outlier Detection ─────────────────────────────────────────
def detect_outliers_iqr(data: np.ndarray, k: float = 1.5) -> np.ndarray:
    """IQR-based outlier detection. Returns boolean mask (True = outlier)."""
    Q1, Q3 = np.percentile(data, [25, 75])
    IQR = Q3 - Q1
    return (data < Q1 - k*IQR) | (data > Q3 + k*IQR)

def detect_outliers_zscore(data: np.ndarray, threshold: float = 3.0) -> np.ndarray:
    """Z-score outlier detection."""
    z = np.abs(stats.zscore(data))
    return z > threshold

# ── Drift Correction ──────────────────────────────────────────
def polynomial_detrend(signal_values: np.ndarray, degree: int = 2) -> np.ndarray:
    """Remove polynomial trend (drift) from signal."""
    t = np.arange(len(signal_values))
    coeffs = np.polyfit(t, signal_values, degree)
    trend = np.polyval(coeffs, t)
    return signal_values - trend

# ── Synchronization of Multiple Sensors ──────────────────────
def synchronize_sensors(sensor_dfs: list, target_freq_hz: float) -> pd.DataFrame:
    """Resample and align multiple sensor streams to common time axis."""
    period = f"{int(1000/target_freq_hz)}ms"
    resampled = []
    for df in sensor_dfs:
        df.index = pd.to_datetime(df.index)
        resampled.append(df.resample(period).interpolate("linear"))
    return pd.concat(resampled, axis=1).dropna()

6. Finite Element Analysis Concepts

Key Checks for FEA Results

# Mesh quality metrics (compute from mesh data)
def aspect_ratio(element_lengths):
    """Aspect ratio: max_edge / min_edge. Should be < 5 for good mesh."""
    return max(element_lengths) / min(element_lengths)

# Convergence study: run at multiple mesh refinements
mesh_sizes = [0.1, 0.05, 0.025, 0.01]  # element size
max_stress = [245, 312, 341, 347]       # MPa (example)

# Plot convergence: result should plateau
import numpy as np
rel_change = np.abs(np.diff(max_stress)) / np.abs(max_stress[:-1]) * 100
print(f"Change with mesh refinement (%): {rel_change}")
# Accept when change < 5% with further refinement

FEA best practices:

• Always run a mesh convergence study (refine until solution changes < 5%)
• Check reaction forces balance applied loads (global equilibrium check)
• Stress singularities at sharp corners → use fillets or accept as mathematical artifact
• Symmetry boundary conditions reduce model size but must match actual symmetry
• Report: element type, mesh size at critical regions, solver tolerance, convergence criterion