Le critère de Routh-Hurwitz : une pierre angulaire de l'analyse des systèmes de contrôle

Le critère de Routh-Hurwitz est une méthode fondamentale de contrôle technique utilisée pour déterminer la stabilité d'un système basé sur son équation caractéristique. Traditionnellement, effectuer ces contrôles manuellement peut être long et sujet à des erreurs, surtout à mesure que l'ordre du système augmente. Automatiser le processus avec des scripts Python améliore significativement l'efficacité et la précision, permettant aux ingénieurs d'analyser rapidement les polynômes de haut ordre, d'intégrer les contrôles de stabilité dans les boucles de conception automatisées et de réduire le risque d'erreurs de calcul.

Cet article fournit un guide complet pour automatiser les contrôles de stabilité de Routh-Hurwitz en utilisant Python. Il couvre le fond mathématique, les détails d'implémentation, la gestion des cas spéciaux et l'intégration pratique dans les workflows d'ingénierie. D'ici la fin, vous serez équipé pour écrire des scripts robustes qui peuvent traiter plusieurs configurations de système, gérer les coefficients symboliques et produire des évaluations de stabilité claires.

C'est quoi le critère de Routh-Hurwitz ?

Nommé d'après Edward John Routh et Adolf Hurwitz, le critère fournit des conditions nécessaires et suffisantes pour la stabilité d'un système linéaire invariant dans le temps (LTI).

(avec a0 > 0)

La méthode construit le tableau de Routh à partir des coefficients. Le système est stable si, et seulement si, tous les éléments de la première colonne du tableau sont positifs. Tout changement de signe indique des pôles instables, et un zéro dans la première colonne ou une ligne entière de zéros indique une stabilité marginale ou la présence de racines symétriquement localisées.

Le critère est largement utilisé car il évite de résoudre le polynôme lui-même et fonctionne directement avec des coefficients. Cependant, la construction manuelle pour les polynômes de degré 5 ou plus devient fastidieuse et sujette à erreur, rendant l'automatisation très précieuse.

Pourquoi Automatiser les contrôles de stabilité?

L'automatisation de la procédure Routh-Hurwitz offre des avantages substantiels tant dans le milieu universitaire que dans l'industrie:

  • Speed et efficacité:[ Un script peut évaluer des dizaines de polynômes par seconde, permettant une itération rapide pendant le réglage du contrôleur ou l'identification du système.
  • Acquiescement:[ Élimine les erreurs arithmétiques qui se produisent couramment lors de la construction du tableau à la main, surtout lorsqu'il s'agit de coefficients fractionnels ou symboliques.
  • Les ingénieurs peuvent systématiquement varier les gains, les constantes de temps ou d'autres paramètres de conception et voir immédiatement l'effet sur les limites de stabilité.
  • Intégration avec des simulations plus grandes: Le contrôle de stabilité peut être intégré dans des boucles d'optimisation, des analyses Monte Carlo ou des scripts de conception automatisés.
  • Reproductibilité: Les scripts fournissent une analyse documentée, contrôlée par version, qui peut être facilement partagée et vérifiée.

En déchargeant le calcul de routine à Python, les ingénieurs peuvent se concentrer sur les décisions de conception de niveau supérieur et l'interprétation des résultats.

Fondations de la construction de Routh Array

Avant d'écrire le code, il est essentiel de comprendre l'algorithme que le script doit suivre.

  1. Disposer les coefficients en deux lignes : la première ligne contient des coefficients de puissance égale (à partir de la puissance la plus élevée), et la deuxième ligne contient des coefficients de puissance impair. Par exemple, pour un polynôme du cinquième ordre a0s5 + a1s4 + a2s3 + a3s2 + a4s + a5, la première ligne est [a0, a2, a4] et la deuxième ligne est [a1, a3, a5].
  2. Lignes de pad avec zéros pour assurer une longueur égale si nécessaire.
  3. Calculer les lignes suivantes en utilisant la formule:

(où a1 est le premier élément de la ligne précédente)

Répéter jusqu'à ce que le tableau ait des lignes n+1.

Une nuance cruciale : si un zéro apparaît dans la première colonne d'une ligne, la formule standard échoue. La construction du tableau de Routh doit traiter des cas particuliers :

  • Zero dans la première colonne (mais pas tous les zéros dans la ligne): Remplacer le zéro par un petit nombre positif ε, poursuivre la construction, puis examiner les signes comme ε → 0.
  • Entrez la ligne de zéros:[ Cela indique la présence de racines placées symétriquement (p. ex., conjuguées complexes sur l'axe imaginaire, ou paires de racines avec des signes opposés). Le polynôme auxiliaire formé à partir de la ligne au-dessus de la ligne zéro doit être utilisé pour poursuivre le tableau.

Un script d'automatisation robuste doit détecter et traiter de façon appropriée les deux cas.

Mise en œuvre de l'algorithme Routh-Hurwitz en Python

Bibliothèques et configuration

Nous utiliserons SymPy pour les mathématiques symboliques et NumPy pour les calculs numériques. SymPy est particulièrement utile lorsque les coefficients impliquent des paramètres (p. ex. gains inconnus) qui ne sont pas numériques.

Mise en œuvre numérique de base

Le script le plus simple accepte une liste de coefficients numériques (float ou entier) et construit le tableau Routh en utilisant l'arithmétique en virgule flottante. Ci-dessous est une version mise à jour et élargie du code de base:

import numpy as np

def routh_hurwitz_numeric(coeffs):
 """
 Construct the Routh array for a polynomial with numeric coefficients.
 coeffs: list of coefficients from highest power down (a0, a1, ..., an)
 Returns a tuple (array, stability: str) or raises ValueError if first element is zero.
 """
 if coeffs[0] <= 0:
 raise ValueError("Coefficient a0 must be positive for standard Routh-Hurwitz.")
 n = len(coeffs) - 1
 # Build first two rows
 row1 = np.array([coeffs[i] for i in range(0, n+1, 2)], dtype=float)
 row2 = np.array([coeffs[i] for i in range(1, n+1, 2)], dtype=float)
 # Pad to same length
 max_len = max(len(row1), len(row2))
 row1 = np.pad(row1, (0, max_len - len(row1)))
 row2 = np.pad(row2, (0, max_len - len(row2)))
 routh = [row1.tolist(), row2.tolist()]

 for i in range(2, n+1):
 if routh[-1][0] == 0:
 # Special case: zero in first column
 # Replace with small epsilon (here we modify the row)
 routh[-1][0] = 1e-10 # Use a tiny positive number
 # Mark that epsilon was used (for sign analysis)
 # For simplicity, we assume epsilon > 0; later we can refine
 row = []
 for j in range(len(routh[0]) - 1):
 a = routh[i-2][0]
 b = routh[i-2][j+1] if j+1 < len(routh[i-2]) else 0
 c = routh[i-1][0]
 d = routh[i-1][j+1] if j+1 < len(routh[i-1]) else 0
 if c == 0:
 row.append(0) # Not reached because we handled epsilon above
 else:
 row.append((c * b - a * d) / c)
 if all(abs(x) < 1e-12 for x in row):
 # Entire row of zeros -> handle auxiliary polynomial
 return handle_auxiliary_row(routh, row, coeffs)
 routh.append(row)
 # Check stability
 first_col = [routh[i][0] for i in range(n+1)]
 if all(x > 0 for x in first_col):
 return routh, "Stable"
 else:
 return routh, "Unstable"

def handle_auxiliary_row(routh, zero_row, coeffs):
 # Extract row above zero row (the one that generated the auxiliary polynomial)
 prev_row = routh[-1]
 # Build auxiliary polynomial from prev_row: s^? (using degrees)
 # This is a simplified placeholder; full implementation requires polynomial differentiation
 # For detail, see reference.
 # For now, we raise an error indicating advanced handling needed.
 raise NotImplementedError("Entire row of zeros: auxiliary polynomial method required. Consider using sympy implementation.")

Bien que cette implémentation numérique fonctionne pour de nombreux cas, elle devient fragile près des singularités. L'approche de remplacement ε nécessite un suivi attentif des changements de signe.

Mise en œuvre symbolique en utilisant SymPy

SymPy , les capacités rationnelles arithmétiques et limites permettent une manipulation correcte et précise des cas de ligne zéro et auxiliaire. Voici une version symbolique complète:

import sympy as sp

def routh_hurwitz_symbolic(coeffs):
 """
 coeffs: list of symbolic or numeric coefficients (a0, a1, ..., an), a0 > 0.
 Returns the Routh array (list of lists) and a stability message.
 """
 coeffs = [sp.sympify(c) for c in coeffs]
 n = len(coeffs) - 1
 # First two rows
 row1 = [coeffs[i] for i in range(0, n+1, 2)]
 row2 = [coeffs[i] for i in range(1, n+1, 2)]
 # Pad
 max_len = max(len(row1), len(row2))
 row1 += [0] * (max_len - len(row1))
 row2 += [0] * (max_len - len(row2))
 routh = [row1, row2]

 epsilon = sp.symbols('epsilon', positive=True)

 for i in range(2, n+1):
 prev_row1 = routh[i-2]
 prev_row2 = routh[i-1]
 first_prev = prev_row2[0]
 # Check for zero in first column
 if first_prev == 0:
 if all(x == 0 for x in prev_row2):
 # Entire row of zeros
 # Build auxiliary polynomial from previous row
 # Ex: if row above zero row is [a, b, c, ...] -> auxiliary poly: a*s^? + b*s^? + ...
 # Need to reconstruct degrees. Function robuster needed.
 # For brevity, we refer to the extended implementation.
 return routh, "Marginal stability detected; auxiliary polynomial needed"
 else:
 # Replace zero with epsilon
 prev_row2 = [epsilon if j == 0 else prev_row2[j] for j in range(len(prev_row2))]
 first_prev = epsilon
 new_row = []
 for j in range(len(prev_row1) - 1):
 a = prev_row1[0]
 b = prev_row1[j+1] if j+1 < len(prev_row1) else 0
 c = first_prev
 d = prev_row2[j+1] if j+1 < len(prev_row2) else 0
 if c == 0:
 value = 0
 else:
 value = (c * b - a * d) / c
 new_row.append(sp.simplify(value))
 # Simplify the row
 new_row = [sp.simplify(x) for x in new_row]
 routh.append(new_row)
 # After building row, if epsilon was used, take limit epsilon -> 0+
 # This simplifies the row to a numeric result if possible.
 if epsilon in sp.flatten([sp.preorder_traversal(x) for x in new_row]):
 new_row = [sp.limit(x, epsilon, 0) for x in new_row]
 routh[-1] = new_row

 # Check first column signs
 first_col = [routh[i][0] for i in range(n+1)]
 # If any symbol still present, cannot decide numerically; user must substitute.
 if any(sp.sympify(x).has(sp.Symbol) for x in first_col):
 return routh, "Indeterminate due to symbolic parameters; substitute numeric values."
 signs = [sp.sign(x) for x in first_col]
 if all(s == 1 for s in signs):
 return routh, "Stable"
 elif any(s == -1 for s in signs):
 return routh, "Unstable"
 else:
 return routh, "Marginal stability (zeros in first column)"

Cette version symbolique gère correctement les zéros et, avec une logique supplémentaire, peut gérer des lignes entières de zéros. Elle est idéale pour les analyses paramétrées où les coefficients contiennent des variables symboliques comme .

Manipulation de rangées entières de zéros

Lorsqu'une rangée de zéros apparaît, la construction standard doit passer à la méthode polynôme auxiliaire. Le polynôme auxiliaire est formé à partir de la rangée directement au-dessus de la rangée zéro. Les coefficients de cette rangée correspondent à des puissances égales de s si la rangée zéro est à un indice impair, etc. Les dérivés du polynôme auxiliaire fournissent la ligne de remplacement. La mise en œuvre complète de la manipulation nécessite plusieurs lignes supplémentaires; le module SymPy polys[ peut aider. Pour des raisons d'espace, nous référons l'exemple complet dans le code source SymPy ou la Python Control Systems Library[ qui comprend les fonctions de Routh-Hurwitz intégrées.

Essais et validation

Tout script automatisé doit être testé en fonction de cas connus.

  • Polynômes stables (p. ex. -> stables)
  • Polynômes instables (p. ex. -> instables en raison du changement de signe)
  • Polynômes avec zéro dans la première colonne (p. ex. -> marginal ou instable selon les racines)
  • Polynômes avec une rangée entière de zéros (p. ex. -> marginal)
  • Polynômes de haut ordre (p. ex. degré 10) pour vérifier les performances.

Comparer les résultats avec le calcul manuel ou la sortie connue de la bibliothèque de contrôle ou des fonctions Routh dédiées.

Intégration dans les flux de travail en génie

Une fois le script fiable, l'intégrer dans un environnement Python plus large :

Balayage et plis de paramètres

Utilisez NumPy ou pandas pour générer une grille de valeurs de paramètres (p. ex., gagner K de 0 à 100). Pour chaque valeur, calculez les coefficients polynômes caractéristiques (via la fonction de transfert système), exécutez le contrôle de stabilité et stockez le résultat. Visualisez les régions de stabilité avec matplotlib:

import numpy as np
import matplotlib.pyplot as plt
from control import tf, feedback

def check_stability_for_gain(K):
 # Example: unity feedback with plant G(s) = K/(s^3 + 3s^2 + 2s)
 G = tf([K], [1, 3, 2, 0])
 T = feedback(G, 1)
 poly = T.den[0][0] # denominator coefficients
 stable = routh_hurwitz_numeric(poly) # call your function
 return stable

Ks = np.linspace(0, 20, 100)
stable_list = [check_stability_for_gain(K) for K in Ks]
plt.plot(Ks, stable_list)
plt.xlabel('Gain K')
plt.ylabel('Stable (1) / Unstable (0)')
plt.show()

Ces placettes révèlent rapidement des marges de stabilité (p. ex., la marge de gain où la stabilité change).

Optimisation de la conception automatisée

Intégrer le contrôle de stabilité comme une contrainte dans l'optimisation. Par exemple, utilisez scipy.optimize pour minimiser une fonction de coût tout en exigeant le critère Routh-Hurwitz pour retourner "Stable". La version symbolique permet une évaluation sans gradient.

Intégration avec les carnets Jupyter

Les cahiers Jupyter sont parfaits pour l'analyse interactive. Combinez la fonction Routh-Hurwitz avec la jolie impression de Sympy pour afficher le tableau étape par étape, ce qui le rend éducatif et débogable.

Considérations avancées

Précision numérique

Pour les coefficients flottants, l'algorithme peut souffrir d'erreurs d'annulation. Utilisez SymPy avec des nombres rationnels lorsque c'est possible, ou utilisez des flotteurs de haute précision via . Sinon, le module numpy.polynomial peut calculer les racines directement, ce qui est plus simple pour les polynômes numériques.

Manipulation des paramètres symboliques avec hypothèses

Lorsque les coefficients impliquent des symboles (p. ex., K, -), le script peut produire des expressions nécessitant une analyse manuelle des signes. Utilisez SymPy= avec des hypothèses pour simplifier en fonction de gammes connues (p. ex., K > 0). Cela peut partiellement automatiser la décision de stabilité.

Parallélisation

Pour les grands balayages de paramètres, parallélisez en utilisant ou . Chaque contrôle de stabilité est indépendant.

Conclusion et pratiques exemplaires

Automatiser les contrôles de stabilité de Routh-Hurwitz avec les scripts Python transforme un processus manuel fastidieux en un outil rapide, fiable et extensible.

  • Commencez par une compréhension claire de l'algorithme, y compris les cas spéciaux.
  • Utiliser SymPy[ pour les coefficients symboliques et la manipulation exacte du cas de ligne zéro.
  • Utiliser NumPy pour les polynômes numériques purs lorsque la vitesse est primordiale.
  • Testez votre implémentation avec une variété de polynômes.
  • Intégrer la fonction dans des pipelines d'analyse plus larges pour la conception et l'optimisation.

En suivant les exemples et les lignes directrices de cet article, vous pouvez créer des scripts d'automatisation robustes qui servent de pierre angulaire à votre travail de systèmes de contrôle. Le temps économisé vous permettra d'explorer plus d'alternatives de conception et d'obtenir de meilleures performances du système.

Note:[ Le code prévu dans le présent article est à des fins éducatives. Pour l'utilisation de la production, envisager de contribuer à des paquets établis comme Python Control Systems Library, qui comprend des fonctions d'analyse de stabilité bien testées de Routh-Hurwitz et d'autres.

Avec les scripts en main, vous êtes maintenant prêt à automatiser les contrôles de stabilité pour les systèmes de n'importe quel ordre, libérant votre énergie mentale pour les défis créatifs de la conception de système de contrôle.