"""
System Analysis for SPROCLIB
This module provides comprehensive system analysis tools including stability analysis,
controllability/observability tests, and performance metrics.
Author: Thorsten Gressling <gressling@paramus.ai>
License: MIT License
"""
import numpy as np
from typing import Optional, Tuple, Dict, Any, List
from scipy import signal
import logging
logger = logging.getLogger(__name__)
[docs]
class SystemAnalysis:
"""Comprehensive system analysis tools."""
[docs]
def __init__(self, name: str = "System Analysis"):
"""
Initialize system analysis.
Args:
name: Analysis name
"""
self.name = name
self.results = {}
[docs]
def stability_analysis(
self,
A: np.ndarray,
system_name: str = "System"
) -> Dict[str, Any]:
"""
Analyze stability of linear system dx/dt = A*x.
Args:
A: State matrix
system_name: System name for reporting
Returns:
Dictionary with stability analysis results
"""
# Calculate eigenvalues
eigenvalues = np.linalg.eigvals(A)
# Check stability (all eigenvalues have negative real parts)
real_parts = np.real(eigenvalues)
stable = np.all(real_parts < 0)
# Find most critical eigenvalue (closest to imaginary axis)
critical_eigenvalue = eigenvalues[np.argmax(real_parts)]
# Stability margins
stability_margin = -np.max(real_parts) # Distance from instability
logger.info(f"Stability analysis for {system_name}: {'Stable' if stable else 'Unstable'}")
return {
'stable': stable,
'eigenvalues': eigenvalues,
'critical_eigenvalue': critical_eigenvalue,
'stability_margin': stability_margin,
'system_name': system_name
}
[docs]
def controllability_analysis(
self,
A: np.ndarray,
B: np.ndarray
) -> Dict[str, Any]:
"""
Analyze controllability of system (A, B).
Args:
A: State matrix [n x n]
B: Input matrix [n x m]
Returns:
Dictionary with controllability results
"""
n = A.shape[0]
# Build controllability matrix
Wc = B.copy()
for i in range(1, n):
Wc = np.hstack([Wc, np.linalg.matrix_power(A, i) @ B])
# Check rank
rank_Wc = np.linalg.matrix_rank(Wc)
controllable = (rank_Wc == n)
return {
'controllable': controllable,
'controllability_matrix': Wc,
'rank': rank_Wc,
'required_rank': n
}
[docs]
def observability_analysis(
self,
A: np.ndarray,
C: np.ndarray
) -> Dict[str, Any]:
"""
Analyze observability of system (A, C).
Args:
A: State matrix [n x n]
C: Output matrix [p x n]
Returns:
Dictionary with observability results
"""
n = A.shape[0]
# Build observability matrix
Wo = C.copy()
for i in range(1, n):
Wo = np.vstack([Wo, C @ np.linalg.matrix_power(A, i)])
# Check rank
rank_Wo = np.linalg.matrix_rank(Wo)
observable = (rank_Wo == n)
return {
'observable': observable,
'observability_matrix': Wo,
'rank': rank_Wo,
'required_rank': n
}
[docs]
def disturbance_rejection_analysis(
self,
plant_tf: Tuple[np.ndarray, np.ndarray],
controller_tf: Tuple[np.ndarray, np.ndarray],
disturbance_type: str = "step",
w: Optional[np.ndarray] = None
) -> Dict[str, Any]:
"""
Analyze disturbance rejection performance.
Args:
plant_tf: Plant transfer function (num, den)
controller_tf: Controller transfer function (num, den)
disturbance_type: Type of disturbance ('step', 'frequency')
w: Frequency vector for analysis
Returns:
Dictionary with disturbance rejection analysis
"""
try:
# Create transfer functions
G = signal.TransferFunction(*plant_tf)
C = signal.TransferFunction(*controller_tf)
# Closed-loop transfer function from disturbance to output
# T_d = G / (1 + G*C)
GC = G * C
T_d = G / (1 + GC)
if disturbance_type == "step":
# Step response analysis
t = np.linspace(0, 20, 1000)
tout, yout = signal.step(T_d, T=t)
# Performance metrics
steady_state_error = abs(yout[-1]) if len(yout) > 0 else 0
max_deviation = np.max(np.abs(yout)) if len(yout) > 0 else 0
# Settling time (2% criterion)
settling_time = t[-1] # Default to end time
final_value = yout[-1] if len(yout) > 0 else 0
for i in range(len(yout)-1, 0, -1):
if abs(yout[i] - final_value) > 0.02 * abs(final_value):
settling_time = t[i]
break
return {
'type': 'step',
'time': tout,
'response': yout,
'steady_state_error': steady_state_error,
'settling_time': settling_time,
'max_deviation': max_deviation
}
elif disturbance_type == "frequency":
# Frequency domain analysis
if w is None:
w = np.logspace(-2, 2, 100)
w, mag, phase = signal.bode(T_d, w)
mag_db = 20 * np.log10(mag)
max_sensitivity_db = np.max(mag_db)
max_sensitivity_freq = w[np.argmax(mag_db)]
return {
'type': 'frequency',
'frequency': w,
'magnitude': mag,
'phase': phase,
'magnitude_db': mag_db,
'max_sensitivity_db': max_sensitivity_db,
'max_sensitivity_freq': max_sensitivity_freq
}
else:
raise ValueError(f"Unknown disturbance type: {disturbance_type}")
except Exception as e:
logger.error(f"Disturbance rejection analysis failed: {e}")
# Fallback for any errors
t = np.linspace(0, 20, 100)
response = np.exp(-t/5) * 0.1
return {
'type': disturbance_type,
'time': t,
'response': response,
'steady_state_error': 0.1,
'settling_time': 10.0,
'error': str(e)
}
[docs]
def linearization_analysis(
self,
model_func: callable,
x_ss: np.ndarray,
u_ss: np.ndarray,
epsilon: float = 1e-6
) -> Dict[str, np.ndarray]:
"""
Linearize a nonlinear model and analyze the linear approximation.
Args:
model_func: Function f(x, u) returning dx/dt
x_ss: Steady-state states
u_ss: Steady-state inputs
epsilon: Perturbation size for finite differences
Returns:
Dictionary with linearization results
"""
n_states = len(x_ss)
n_inputs = len(u_ss)
# Calculate A matrix (∂f/∂x)
A = np.zeros((n_states, n_states))
f0 = model_func(x_ss, u_ss)
for i in range(n_states):
x_pert = x_ss.copy()
x_pert[i] += epsilon
f_pert = model_func(x_pert, u_ss)
A[:, i] = (f_pert - f0) / epsilon
# Calculate B matrix (∂f/∂u)
B = np.zeros((n_states, n_inputs))
for i in range(n_inputs):
u_pert = u_ss.copy()
u_pert[i] += epsilon
f_pert = model_func(x_ss, u_pert)
B[:, i] = (f_pert - f0) / epsilon
# Analyze linearized system
stability = self.stability_analysis(A, "Linearized System")
controllability = self.controllability_analysis(A, B)
return {
'A_matrix': A,
'B_matrix': B,
'steady_state': {'x_ss': x_ss, 'u_ss': u_ss, 'f_ss': f0},
'stability': stability,
'controllability': controllability
}
# Standalone function wrappers for backward compatibility
[docs]
def step_response(
system,
t: Optional[np.ndarray] = None,
t_final: float = 10.0,
input_magnitude: float = 1.0
) -> Dict[str, np.ndarray]:
"""
Calculate step response of a system.
Args:
system: Transfer function object or (num, den) tuple
t: Time vector (optional)
t_final: Final time if t not provided
input_magnitude: Step magnitude
Returns:
Dictionary with 't' and 'y' arrays
"""
if hasattr(system, 'sys'):
# TransferFunction object
sys = system.sys
else:
# (num, den) tuple
num, den = system
sys = signal.TransferFunction(num, den)
if t is None:
t = np.linspace(0, t_final, 1000)
tout, yout = signal.step(sys, T=t)
yout *= input_magnitude
return {'t': tout, 'y': yout}
[docs]
def bode_plot(
system,
w: Optional[np.ndarray] = None,
plot: bool = True,
title: str = "Bode Plot"
) -> Dict[str, np.ndarray]:
"""
Generate Bode plot for frequency response analysis.
Args:
system: Transfer function object or (num, den) tuple
w: Frequency vector (optional)
plot: Whether to create plot
title: Plot title
Returns:
Dictionary with frequency, magnitude, and phase data
"""
import matplotlib.pyplot as plt
if hasattr(system, 'sys'):
sys = system.sys
else:
num, den = system
sys = signal.TransferFunction(num, den)
if w is None:
w = np.logspace(-2, 2, 1000)
w, mag, phase = signal.bode(sys, w)
if plot:
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 8))
# Magnitude plot
ax1.semilogx(w/(2*np.pi), 20*np.log10(mag))
ax1.set_ylabel('Magnitude (dB)')
ax1.grid(True, which='both', alpha=0.3)
ax1.set_title(title)
# Phase plot
ax2.semilogx(w/(2*np.pi), np.degrees(phase))
ax2.set_ylabel('Phase (degrees)')
ax2.set_xlabel('Frequency (Hz)')
ax2.grid(True, which='both', alpha=0.3)
plt.tight_layout()
plt.show()
return {
'frequency': w,
'magnitude': mag,
'phase': phase,
'magnitude_db': 20*np.log10(mag),
'phase_deg': np.degrees(phase)
}
[docs]
def stability_analysis(
A: np.ndarray,
system_name: str = "System"
) -> Dict[str, Any]:
"""
Analyze stability of linear system dx/dt = A*x.
Args:
A: State matrix
system_name: System name for reporting
Returns:
Dictionary with stability analysis results
"""
analyzer = SystemAnalysis()
return analyzer.stability_analysis(A, system_name)