Sciences Industrielles pour l’Ingénieur(e)

Cycle 3 : Méthodes Numériques

DM : Outils numériques pour la

résolution d’équations différentielles non
linéaires – Corrigé
 Méthodes Numériques

Introduction
Ce DM propose plusieurs exercices qui font le lien entre des lois de comportement établies analytiquement
en classe et les méthodes numériques associées permettant de tracer les courbes issues de ces lois.

Les méthodes évoquées dans ces exercices sont à connaître. Toutefois il est hors de question d’apprendre le
code par cœur, c’est l’idée derrière qui est à retenir.

Exercice 1 : Tracer une loi de comportement issue de la dynamique par le

schéma d’Euler explicite
I Problème : Centrifugeuse de laboratoire
I.1 Mise en situation
On considère une centrifugeuse d’éprouvette de laboratoire. Elle est constituée du bâti 0, d’un bras 1 et d’une
éprouvette 2. L’éprouvette contient deux liquides de masses volumiques différentes. Par rotation rapide du
bras, l’éprouvette s’incline, et la force centrifuge sépare les deux liquides (le liquide de masse volumique la
plus lourde étant projeté au fond de l’éprouvette).

(a) Centrifugeuse (b) Manipulation des éprouvettes

Figure 1 – Centrifugeuse de laboratoire

I.2 Cahier des charges

« requirement » « requirement »
« requirement »
Vitesse angulaire Sécurité
Centrifugeuse
Id="1.2" Id="1.1"
Id="1"
Text="La vitesse angulaire doit être Text="Le fonctionnement de la cen-
Text="Séparer deux liquides."
suffisante." trifugeuse doit être sécure."

« requirement »
« requirement »
« requirement » Dimensions
Accélération
Stabilité
Id="1.4"
Id="1.3"
Id="1.5" Text="Les dimensions de la centri-
Text="L’accélération minimale du
Text="La centrifugeuse doit être fugeuse doivent être adaptées à une
centre de gravité de l’éprouvette
stable." implantation en laboratoire."
doit être d’au moins 24 m/s2."

Figure 2 – Extrait du diagramme des exigences de la centrifugeuse

DM PSI* 1/34
 Méthodes Numériques

I.3 Paramétrage
On définit :
— le repère R0 = (O, − →
x0 , −
→
y0 , −
→
z0 ) est lié au bâti 0 ;
!
z0
— le repère R1 = (O, x1 , y1 , −
−→ −
→ →
z1 ) est lié au bras 1 de masse né-
gligeable. Les solides 0 et 1 sont en liaison pivot d’axe (O, − →
z0 ) 1
et
(−
→x0 , −
→
x1 ) = ( −
→
y0 , −
→
y1 ) = α ;
O⇥ A
— le repère R = ( A, −
2
→
x ,−
2
→
y ,−
2
→
z ) est lié à l’éprouvette 2 assimilée
2
à une masse ponctuelle m = 20g en G. Les solides 1 et 2 sont
en liaison pivot d’axe ( A, − →
y1 ) et
−→ −
→ −→ −
→
( x1 , x2 ) = ( z1 , z2 ) = β. 2 ⇥
−→ −→ 0 G
De plus, OA = a− →
x1 et GA = −b− →
x2 avec a = 10cm et b = 13cm.
L’accélération de la pesanteur sera notée − →g = − g− →
z avec 0
g = 9, 81m/s². Enfin la liaison pivot entre 0 est 1 est motorisée,
le couple fourni à 1 par le moteur est Cm −→
z0 . En outre la liaison
pivot en A est sujette à du frottement visqueux de constante
µ = 4 · 10−3 N.m.s/rad.

II Étude dynamique
Question 1 – Montrer que les équations du mouvement qui régissent le système s’écrivent :
 
m α̈ ( a + b cos( β)) − 2bα̇ β̇ sin( β) · ( a + b cos( β)) = Cm (1)

mb2 β̈ + µ β̇ + mb ( a + b cos( β)) α̇2 sin( β) = mgb cos( β) (2)

Correction On commence par tracer le graphe des liaisons avec les AM autres que celles transmises par
les liaisons et on dessine les figures de changement de base pour les calculs vectoriels :

−mg−
→
z0
Pivot (O, −
→
z0 ) Pivot ( A, −
→
y1 )
0 1 2

Cm −
→
z0 frottement visqueux −µ β̇−
→
y1

−
→ −
→ −
→ −
→
y1 y0 x2 x1
α β
−
→ −
→
z2
x1
α −
→ β −
→
x0 −
→ z1
−
→
z0 = −
→
z1 y1 = −
→
y2

On a clairement une chaîne ouverte, la stratégie à adopter est donc immédiate. On procède alors de l’iso-
lement impliquant le moins de calculs à l’isolement impliquant le plus de calculs.

Remarque : même pour une chaîne ouverte aussi simple que celle-ci, on va voir qu’il faut déjà être très
attentif aux étapes de calculs pour ne pas se tromper.

DM PSI* 2/34
 Méthodes Numériques

Isolement de 2 :

On isole {2}, le Bilan des Actions Mécaniques Extérieures (BAME) est déjà réalisé sur le graphe des liaisons.
On profite du zéro de liaison pivot et on applique le Théorème du Moment Dynamique en A projeté selon
−
→
y1 = −
→y2 (vecteur principal du pivot en A) :

−−−−−→ − −→
δ( A, 2/0) · →
y1 = AG ∧ (−mg)−
→
z0 · −
→
y1 − µ β̇−
→
y1 · −
→

y1 = mgb cos( β) − µ β̇

−−−−−→ →
Calcul de δ( A, 2/0) · −
y1 :
−−−−−→ − −−−−−→ −→ −−−−−−→ →
δ( A, 2/0) · →
y1 = δ( G, 2/0) + AG ∧ mΓ( G, 2/0) · −
y1
Varignon
−−−−−−→
−−−−−→ dσ ( G, 2/0) −−−−−−→ − →
Or δ( G, 2/0) = car G est centre d’inertie de 2 et σ ( G, 2/0) = 0 (2 est masse ponctuelle en G).
dt
0
Enfin on calcule :
−→ h−−−−→
dOG
= aα̇−
→
y1 + b Ω(2/0) ∧ −
→
x2 = α̇ [ a + b cos( β)] −
→
y1 − b β̇−
→
i
z2
dt
0
−−−−−−→ 
Γ( G, 2/0) = α̈ ( a + b cos( β)) − bα̇ β̇ sin( β) − bα̇ β̇ sin( β) −
→
y1 − α̇2 [ a + b cos( β)] −
→
x1 − b β̇2 −
→
x2 − b β̈−
→
z2
Donc
−→ on−−−−−−
déduit→
AG ∧ mΓ( G, 2/0) = mb α̈ ( a + b cos( β)) − 2bα̇ β̇ sin( β) −
→
z2 + mbα̇2 [ a + b cos( β)] sin( β)−
→
y1 + mb2 β̈−
→

y1 .

−−−−−→ → −→ −−−−−−→ →

Enfin δ( A, 2/0) · −
y1 = AG ∧ mΓ( G, 2/0) · −
y1 = mbα̇2 [ a + b cos( β)] sin( β) + mb2 β̈ .

On retrouve l’équation (2).

Isolement de {1+2} :

On isole {1+2}, le BAME est déjà réalisé sur le graphe des liaisons. On profite du zéro de liaison pivot et
on applique le Théorème du Moment Dynamique en O projeté selon − →
z0 = −→
z1 (vecteur principal du pivot
en O) :

−−−−−→ − −−−−−→ → −→

δ(O, 2/0) · →
z0 + δ(O, 1/0) · −
z0 = OG ∧ (−mg)−
→
z0 · −
→

z0 + Cm = Cm

−−−−−→ →
Calcul de δ(O, 1/0) · −
z0 :
−−−−−→ −
δ(O, 1/0) · →
z0 = 0 (masse de 1 négligeable).
−−−−−→ →
Calcul de δ(O, 2/0) · −
z0 :

−−−−−→ − −−−−−→ −→ −−−−−−→ →

δ(O, 2/0) · →
z0 = δ( G, 2/0) + OG ∧ mΓ( G, 2/0) · −
z0
Varignon
−−−−−−→
−−−−−→ dσ ( G, 2/0) −−−−−−→ − →
Or δ( G, 2/0) = car G est centre d’inertie de 2 et σ ( G, 2/0) = 0 (2 est masse ponctuelle en G).
dt
0 −→ −−−−−−→ →
ASTUCE : On profite du fait que OG ∧ mΓ( G, 2/0) · − z0 est aussi appelé "produit mixte" et est invariant
par permutation circulaire. Par conséquent :
−→ −−−−−−→ →  −−−−−−→ −  −→
OG ∧ mΓ( G, 2/0) · − z0 = mΓ( G, 2/0) ∧ → z0 · OG

DM PSI* 3/34
 Méthodes Numériques

Remarque : on applique cette astuce car on va pouvoir rentrer le produit vectoriel par −
→
z0 dans la dérivée
−→
seconde du vecteur position OG...
−→ →
d2 OG ∧ −
z0
= α̈ ( a + b cos( β)) − 2bα̇ β̇ sin( β) −
→ h 2
x1 + α̇ ( a + b cos( β)) + b β̈ sin( β) + β̇2 cos( β) − →
  i
Donc y1 .
dt2
0
−−−−−→ →  −−−−−−→ −  −→
Enfin δ(O, 2/0) · −
z0 = mΓ( G, 2/0) ∧ →
 
z0 · OG = m −2bα̇ β̇ sin( β) + α̈ ( a + b cos( β)) · ( a + b cos( β)) .

On retrouve l’équation (1).

III On suppose que le bras 1 atteint le régime permanent très vite

On suppose que le régime permanent est atteint très vite pour la vitesse de rotation α̇. Cela permet alors de
poser les hypothèses suivantes :
• on a α̈ = 0rad/s² et α̇ = 100tr/min pendant tout l’essai, ce qui permet de n’utiliser que l’équation (2) ;
• au départ on a β = 90◦ et β̇ = 0rad/s.
On cherche à tracer l’évolution de β en fonction du temps.

III.1 Avec la fonction odeint de la librairie scipy

On importe la fonction odeint de la librairie scipy :
1 from scipy . integrate import odeint

Dans le Shell, la commande help(odeint) indique que cette fonction résout une équation différentielle de
la forme Ẏ = f (Y, t) où Y peut être un vecteur. Autrement dit cette fonction Python a besoin de connaître la
forme de f ainsi que les conditions initiales.
 
β
On pose Y = . On cherche une fonction f faisant intervenir les composantes du vecteur Y ainsi que le
β̇
temps t.

Question 2 – Exprimer la fonction f (Y, t) à partir de l’équation du mouvement (2). Écrire alors une fonc-
tion f(Y:array,t:float)->array qui renvoie le second membre de l’équation Ẏ = f (Y, t). On pourra
renvoyer le résultat sous forme de vecteur [Link].
 
β̇
Correction On a Ẏ = . En utilisant les équations de la dynamique et les hypothèses de l’énoncé on
β̈
obtient :
 
β̇
Ẏ =  µ a + b cos( β) 2 g  = f (Y, t)
− 2 β̇ − α̇ sin( β) − cos( β)
mb b b
En utilisant β = Y [0] et β̇ = Y [1] on peut écrire le script Python suivant :
1 # Importation des librairies
2 import numpy as np
3 from math import cos ,sin ,pi ,sqrt , floor
4 import matplotlib . pyplot as plt
5 from scipy . optimize import fsolve
6 from scipy . integrate import odeint
7
8 plt. close (’all ’) # On ferme les figures déjà ouvertes
9
10 # Données géométriques (en m)
11 a = 0.1
12 b = 0.13

DM PSI* 4/34
 Méthodes Numériques

13
14 # Masse (en kg)
15 m = 0.02
16 # Amortissement (USI)
17 mu = 4e−3
18
19 # Accélération de la pesanteur
20 g = 9.81
21
22 #——————————————————————————-
23 # Avec scipy
24
25 # vitesse angulaire dot(alpha) constante
26 dalpha = 100 * pi /30
27
28 # Question 2
29 def f(Y,t):
30 r e t u r n np. array ([Y[1],−1/b * ( a+b * cos(Y[0])) * dalpha * * 2 * sin(Y[0])−mu /(m * b
* * 2 ) * Y[1]+g/b * cos(Y[0]) ])
31
32 # Conditions initiales sur beta, attention les fonctions trigo prennent des rad
33 CI = np. array ([ pi /2 ,0])

Remarque : le temps t n’intervient pas hors de β ou β̇ pour notre étude, mais il pourrait d’où la méthode
numérique Ẏ = f (Y, t).

Question 3 – En utilisant help(odeint), résoudre l’équation différentielle et tracer la courbe d’évolution

de β ainsi obtenue. On fera l’essai sur 1 seconde avec 1000 points grâce à la commande
temps = [Link](0,1,1000). Montrer que la courbe obtenue est la suivante :

90
80
70
60
angle beta

50
40
30
20
10
0.0 0.2 0.4 0.6 0.8 1.0
temps (s)

Figure 3 – Résultat fourni par scipy pour résoudre l’équation de la dynamique

Correction On écrit le script à la suite du script écrit en question précédente :

1 # Temps
2 N = 1000 # Nombre de points
3 temps = np. linspace (0,1,N)
4

DM PSI* 5/34
 Méthodes Numériques

5 # Calcul de chaque valeur de l’angle beta

6 beta_scipy = odeint (f,CI , temps )
7
8 fig1 = plt. figure (1)
9 # Première coordonnée de la matrice beta_scipy ramenée en degrés
10 [Link](temps , beta_scipy [: ,0] * 180/ pi)
11 plt. xlabel (" temps (s)")
12 plt. ylabel (" angle beta")
13 [Link](’both ’)
14 [Link] ()

On retrouve la figure de l’énoncé. On vérifie évidemment qu’elle est physiquement cohérente (condi-
tions initiales, convergence vers une position de régime permanent cohérente vis-à-vis du paramétrage de
l’énoncé).

III.2 Avec un schéma du type Ẏ = f (Y, t) codé manuellement

Question 4 – Avec le schéma d’Euler explicite pour un pas de temps h, donner la relation de récurrence liant
Yi+1 à Yi et f (Yi , ti ). Yi représente la valeur de Y au temps ti .

Yi+1 − Yi
Correction Le schéma d’Euler explicite s’écrit : Ẏi = . On a alors l’équation de récurrence :
h
Yi+1 = Yi + h · f (Yi , ti )

Question 5 – Écrire une fonction Euler(func,Y0:array,start:float,stop:float,n:int)->array

qui prend en entrées :
• une fonction func ;
• un jeu de conditions initiales Y0 au format [Link] ;
• 2 bornes temporelles start et stop avec start<stop ;
• un entier n qui découpe l’intervalle [start,stop] en n segments (c’est le pas d’Euler pour le calcul
de la dérivée).
Cette fonction renvoie un tableau correspondant à la valeur de Y au temps t=stop.

Correction

1 def Euler (func ,Y0 ,start ,stop ,n):

2 h = (stop − start ) / n # Pas de temps déduit
3 t = start # Initialisation du temps
4 Y = Y0 # Conditions initiales
5 f o r i in range (n) :
6 Y = Y + func(Y, t) * h # Mise à jour de Y
7 t = t + h # Mise à jour de t
8 return Y

Question 6 – Construire le tableau beta_euler qui calcule la valeur de Y pour chaque élément t du tableau
temps. On prendra un pas d’Euler n=1000, start=0 et stop=t.

Correction On écrit le script suivant :

1 n = 1000
2 beta_euler = np. array ([ Euler (f, CI , 0, t, n) f o r t in temps ])

DM PSI* 6/34
 Méthodes Numériques

Remarque Ici le pas de temps h (lié au pas d’Euler n) de la question 4 ne correspond pas au pas de temps
dans le tableau temps (qui vaut 1/1000). En effet à chaque item construit du tableau beta_euler, h est
différent puisque l’intervalle [start,stop] change à chaque fois (il est de plus en plus long).

Question 7 – Tracer sur un même graphe l’angle β obtenu avec scipy et la méthode d’Euler explicite. Les
résultats sont-ils cohérents ?
Correction On écrit le script suivant :

1 fig2 = plt. figure (2)

2 [Link](temps , beta_scipy [: ,0] * 180/ pi , label =’scipy ’)
3 [Link](temps , beta_euler [: ,0] * 180/ pi , label =’Euler ’)
4 plt. xlabel (" temps (s)")
5 plt. ylabel (" angle beta")
6 plt. legend ()
7 [Link](’both ’)
8 [Link] ()

On retrouve la même courbe que dans l’énoncé :

90 scipy
Euler
80
70
60
angle beta

50
40
30
20
10
0.0 0.2 0.4 0.6 0.8 1.0
temps (s)

Question 8 – Refaire le tracé mais pour α̇ = 1000tr/min. Que remarque-t-on ?

Correction On relance le script avec α̇ = 1000tr/min, on obtient la courbe suivante :

DM PSI* 7/34
 Méthodes Numériques

scipy
80 Euler
60
40
angle beta 20
0
20
40
60

0.0 0.2 0.4 0.6 0.8 1.0

temps (s)

La réponse scipy montre des oscillations plus marquées, mais les CI sont les mêmes et on converge à
nouveau vers une position d’équilibre en régime permanent. Avec le schéma d’Euler explicite la réponse
diverge.

Question 9 – Jouer avec le pas d’Euler n pour que les courbes se superposent à nouveau. Que peut-on déduire
du schéma d’Euler explicite pour résoudre une équation sous le format Ẏ = f (Y, t) ?

Correction On augmente la valeur du pas d’Euler en prenant n=10000, on obtient :

scipy
80 Euler
60
40
angle beta

20
0
20
40
60

0.0 0.2 0.4 0.6 0.8 1.0

temps (s)

Le schéma d’Euler explicite est dit conditionnellement stable, autrement dit en fonction des paramètres
du système et du pas d’Euler choisi la réponse peut diverger. Ce n’est pas physiquement cohérent puisque
le frottement visqueux pousse à la stabilité du système.

Une solution est de prendre plus de points de calculs entre 2 pas de temps pour aider à la convergence.

DM PSI* 8/34
 Méthodes Numériques

III.3 En établissant une relation de récurrence directement à partir de l’équation du

mouvement
Cette fois-ci on propose de travailler directement sur l’équation du mouvement. On rappelle que l’approxi-
β − 2β i + β i−1
mation de la dérivée seconde s’écrit : β̈ i ≈ i+1 .
∆t2
Question 10 – En réalisant 2 développements limités à l’ordre 2 bien choisis sur la fonction β(t) et en jouant
sur le pas de temps ∆t, retrouver l’approximation proposée par l’énoncé.

Correction On fait les 2 développements limités suivants :

dβ ∆t2 d2 β
β(t + ∆t) ≈ β(t) + ∆t (t) + (t)
dt 2 dt2

dβ ∆t2 d2 β
β(t − ∆t) ≈ β(t) − ∆t (t) + (t)
dt 2 dt2
Par conséquent en sommant les 2 équations :

d2 β
β(t + ∆t) + β(t − ∆t) ≈ 2β(t) + ∆t2
(t)
dt2
On retrouve l’approximation de la dérivée seconde et c’est un schéma centré.

Question 11 – Montrer alors que l’équation de récurrence sur β s’écrit :

β i+1 = Aβ i − B( a + b cos( β i )) sin( β i ) + C cos( β i ) − Dβ i−1

avec A, B, C et D des paramètres positifs à exprimer et dépendant du pas de temps ∆t.

Correction On discrétise l’équation (2) de la dynamique qui devient :

β i+1 − 2β i + β i−1 β − βi
mb2 + µ i +1 + mb ( a + b cos( β i )) α̇2 sin( β i ) = mgb cos( β i )
∆t2 ∆t
En réarrangeant on obtient β i+1 = Aβ i − B( a + b cos( β i )) sin( β i ) + C cos( β i ) − Dβ i−1 avec :

µ∆t
2+ α̇2 ∆t2 g∆t2 1
A= mb2 ; B=   ; C=   ; D=
µ∆t µ∆t µ∆t µ∆t
1+ b 1+ b 1+ 1+
mb2 mb2 mb2 mb2

On choisit le même pas de temps que dans le tableau temps, soit ∆t = 1/1000.

Question 12 – Tracer l’évolution de l’angle β en fonction du temps. Comparer aux résultats fournis par scipy
et par le schéma d’Euler explicite.

Correction On écrit le script suivant :

1 # Les conditions initiales en position et vitesse se traduisent par

2 beta_rec = np. zeros (N)
3 beta_rec [0] = pi /2
4 beta_rec [1] = pi /2
5
6 dt = 1/N # Pas de temps dans le vecteur ’temps’
7
8 A = (2+ mu * dt /(m * b * * 2 ) )/(1+ mu * dt /(m * b * * 2 ) )
9 B = dalpha * * 2 * dt * * 2 / ( b * (1+ mu * dt /(m * b * * 2 ) ))
10 C = g * dt * * 2 / ( b * (1+ mu * dt /(m * b * * 2 ) ))
11 D = 1/(1+ mu * dt /(m * b * * 2 ) )
12

DM PSI* 9/34
 Méthodes Numériques

13 f o r i in range (2, len( temps )):

14 ajout = A * beta_rec [i−1]−B * ( a+b * cos( beta_rec [i−1])) * sin( beta_rec [i−1]) +
C * cos( beta_rec [i−1]) − D * beta_rec [i−2]
15 beta_rec [i] = ajout
16
17 fig3 = plt. figure (3)
18 [Link](temps , beta_scipy [: ,0] * 180/ pi , label =’scipy ’)
19 [Link](temps , beta_euler [: ,0] * 180/ pi , label =’Euler ’)
20 [Link](temps , beta_rec * 180/ pi , label =’ré currence ’)
21 plt. xlabel (" temps (s)")
22 plt. ylabel (" angle beta")
23 plt. legend ()
24 [Link](’both ’)
25 [Link] ()

On obtient les courbes bien superposées suivantes pour α̇ = 100tr/min et n=10000 pour le pas d’Euler :

90 scipy
Euler
80 récurrence

70
60
angle beta

50
40
30
20
10
0.0 0.2 0.4 0.6 0.8 1.0
temps (s)

DM PSI* 10/34
 Méthodes Numériques

Exercice 2 : Tracer une loi de comportement issue d’un asservissement par

le schéma d’Euler explicite
Remarque On pourra créer un nouveau script Python pour cet exercice pour éviter les notations redon-
dantes avec l’exercice précédent.

I Problème : Élaboration du modèle de connaissance d’un actionneur

linéaire de soutien lombaire placé sur un banc d’essai
Source : d’après Centrale-Supelec MP 2023

I.1 Contexte
On étudie dans ce problème un système de soulagement lombaire utilisant des vérins électriques asservis en
effort. L’orthèse est attachée au patient et les vérins servent à écarter les disques inter-vertébraux et limitent
alors les douleurs dues aux hernies discales par exemple (Figure 1). Un banc d’essai est réalisé pour tester
les actionneurs linéaires (les vérins électriques).

Figure 1 – Modèle numérique biomécanique (à gauche) et détail de la structure vertébrale avec numérotation
des vertèbres (à droite)

I.2 Modèle adopté

Objectif Définir un modèle de connaissance de la dynamique du système permettant d’obtenir les équa-
tions d’un modèle de simulation comparable aux mesures du banc d’essai.

Le système est placé sur un banc d’essai en position horizontale (Figure 2). Dans cette configuration, la
pesanteur est portée par la direction −
→
z 0 . On rappelle que le banc fonctionne selon deux protocoles distincts :
• vérification de la course définie précédemment (non étudiée ici) ;
• vérification de la force exercée par un actionneur linéaire pour effectuer la précontrainte.

DM PSI* 11/34
 Méthodes Numériques

Figure 2 – Banc d’essai

Pour la vérification de la force de précontrainte, l’actionneur linéaire est bloqué en H et un capteur d’effort,
noté cap (Figure 2), supposé indéformable, est placé en K. Cette configuration permet uniquement de valider
la performance relative à la mise en précontrainte.

La raideur du ressort du capteur d’effort de l’actionneur linéaire (Figure 3) a été choisie à partir d’une
campagne d’essais réalisée par différents utilisateurs qui ont exprimé leur ressenti en donnant une note de
confort. L’exploitation des données recueillies a permis au fabricant de déterminer le meilleur compromis
parmi les retours des différents utilisateurs.

Figure 3 – Actionneur lineaire

Dans la configuration spécifique retenue, les points K et H sont immobiles par rapport au châssis du banc
d’essai. Les solides (0) et (4) sont immobiles par rapport au châssis du banc d’essai. L’action mécanique de
l’actionneur linéaire sur le capteur d’effort du banc d’essai est un glisseur de support passant par K et de
−
→
résultante F 4→cap .
−
→
Le capteur du banc d’essai mesure ainsi F 4→cap · −
→
y 0 . On suppose que c’est une image fidèle de la force de
traction exercée par l’actionneur linéaire.

DM PSI* 12/34
 Méthodes Numériques

Figure 4 – Modèle d’étude de l’actionneur linéaire sur le banc d’essai

• On note l (t) le déplacement de (3) par rapport à (4) ;

• on a y(t) = l0 − l (t) avec l0 la longueur à vide du ressort du capteur installé sur le système réel. Le
ressort n’est pas préchargé avant le début de l’essai soit l (t = 0) = l0 ;
• le repère R H, − →x ,−→y ,−
→

0 0 z
0 0lié au châssis du banc d’essai est supposé galiléen.
Les différentes grandeurs utiles à cette partie sont regroupées dans le Tableau 1.

Éléments Caractéristiques et notations

Corps du vérin 0 Masse : m0
Couple moteur : C3→arbre moteur (t) = cm (t)
Moment d’inertie de l’arbre moteur suivant son axe : Im
Moteur
Vitesse de rotation de l’arbre moteur : ωm (t) = ωm/3
Masse négligeable devant les autres masses
Vitesse de rotation en sortie du réducteur : ωr (t) = ω1/3
ωr ( t )
Réducteur planétaire Rapport de réduction : λ =
ωm ( t )
+ pignon 1 Moment d’inertie équivalent reporté sur l’arbre de sortie du réducteur : Ir
Masse négligeable devant les autres masses
Rayon primitif du pignon d’entrée 1 : R1
Transmetteur par engrenage
Rayon primitif du pignon de sortie 2 :
R2 = R1
Moment d’inertie suivant l’axe J, − →
y 0 : IV
Vis + pignon 2 Vis de pas géométrique : pV en m/tr
Masse négligeable devant les autres masses
Raideur : Kres
Ressort
Masse négligeable devant les autres masses
Ensemble 3 Masse : m3 (comprenant la masse du carter réducteur et du carter moteur)
Tige de vérin 4 Masse : m4

Tableau 1 – Caractéristiques principales de l’actionneur linéaire

II Lois de comportement mécaniques

II.1 Transmission du mouvement du moteur vers le ressort
ω1/3
Question 1 – Établir en justifiant la démarche le rapport en fonction de R1 et R2 et faire l’application
ω2/3
numérique.

DM PSI* 13/34
 Méthodes Numériques

Correction Version 1
−−−−−→ − →
Par roulement sans glissement au point I on a : V ( I, 1/2) = 0 .
−−−−−→ −−−−−→
Par composition des vitesses cela donne : V ( I, 1/3) = V ( I, 2/3).
−−−−−→ −−−−−−→ − → −−−−→ − →
Par Varignon : V ( I, 1/3) = V (O, 1/3) + IO ∧ Ω(1/3) = 0 − R1 −
→
x0 ∧ ω1/3 −
→
y0 = − R1 ω1/3 −
→
z0 .
−−−−−→
De la même manière on trouve V ( I, 2/3) = R2 ω2/3 −
→
z0 .

ω1/3 R
Ainsi on a la relation = − 2 = −1 .
ω2/3 R1 AN
Version 2 efficace

ω1/3 R
Les roues dentées tournent autour d’axes fixes dans le repère 3 donc = − 2 = −1 .
ω2/3 R1 AN

−
→
Question 2 – Rappeler le lien entre V J,2/R0 et la vitesse de rotation ω2/R0 −
→
y0 .

Correction Par définition de la liaison hélicoïdale d’axe J, y0

→
: −
− → p
V J,2/R0 = V ω2/R0 −
→
y0.
2π

−
→
Question 3 – Établir en justifiant la démarche (hypothèses et propriétés cinématiques) l’expression de V J,3/R0
en fonction de ωm (t), pV et λ.
−
→ −
→ −
→ −
→ −
→
Correction Par composition des vitesses : V J,3/R0 = V J,3/2 + V J,2/R0 , or V J,3/2 = 0 car J est situé sur
−
→ p
l’axe de rotation de la liaison pivot entre (3) et (2). De plus V J,2/R0 = V ω2/R0 −
→
y 0.
2π
Encore par composition des vitesses ω2/R0 = ω2/3 + ω3/R0 et ω3/R0 = 0 par la liaison glissière entre (3) et
(0). De plus ω2/3 = −ω1/3 = −λωm (t).

−
→ p
Pour conclure V J,3/R0 = − V λωm (t)−
→
y0.
2π

Question 4 – Par composition des vecteurs vitesse en K entre les solides (4), (3) et (0), déterminer la relation
entre ẏ(t), ωm (t), pV et λ.
−
→ −
→ −
→
Correction On a par composition des vitesses V K,4/0 = V K,4/3 + V K,3/0 .
−
→ −
→
Or V K,4/0 = 0 car les solides (0) et (4) sont immobiles.
−
→
De plus V K,4/3 = ẏ(t)−
→
y 0 par définition de y(t).
−
→ −
→ −
→ −→
Enfin V K,3/0 = V K,3/R0 car (0) est immobile, de plus V K,3/R0 = V J,3/R0 par mouvement de translation
−
→ p
de 3/R0 et V J,3/R0 = − V λωm (t)− →y 0 d’après la question précédente.
2π
p
Finalement ẏ(t) = V λωm (t) .
2π

Pour la suite, on définit Ktrans tel que ẏ(t) = Ktrans ωm (t).

II.2 Mesure de l’effort capteur

Question 5 – En prenant soin de préciser le solide isolé et le théorème utilisé, déterminer l’expression littérale
−
→
de la résultante F cap→4 en projection selon −
→
y 0 , en fonction de Kres et y(t).

DM PSI* 14/34
 Méthodes Numériques

Correction On isole le solide (4), il est soumis à :

• le poids du solide (4) selon −
→
z0 ;
• l’action mécanique transmise par la liaison glissière entre 3 et 4 d’axe −
→
y 0;
−→
• l’action du capteur F cap→4 ;
• l’effort de rappel du ressort −K y(t)−→
resy . 0
Le Théorème de la Résultante Statique (solide (4) immobile) projetée selon −
→
y 0 donne :

−
→
0 = F cap→4 · −
→
y 0 − Kres y(t)

II.3 Transmission par engrenages du couple réducteur vers le pignon 2

Le réducteur par engrenages constitué des pignons 1 et 2 possède un rendement unitaire, autrement dit la
puissance en sortie de réducteur (notée P2 ) est exactement celle qui a été fournie en entrée (notée P1 ). On
admet que la puissance mécanique d’un solide en rotation s’écrit sous la forme P = C · ω avec C le couple
qui met en rotation le solide à la vitesse ω.
−−−−→
Question 6 – Établir alors une relation entre le couple en sortie de réducteur Cred→1 = cr (t)−
→
y0 , le couple
−−→ −→ ω1/3
transmis par 1 sur 2 noté C1→2 = c12 (t) y0 ainsi que le rapport de transmission .
ω2/3

Correction Par conservation de la puissance on a P1 = P2 , or P1 = cr (t) · ω1/3 car le pignon 1 est mis en
rotation à la vitesse ω1/3 grâce au couple cr (t) en sortie du réducteur à planétaire.

Aussi P2 = c12 (t) · ω2/3 car le pignon 2 est mis en rotation à la vitesse ω2/3 grâce au couple c12 (t).
ω1/3
Finalement c12 (t) = −cr (t) car = −1.
ω2/3

On admettra de manière similaire que le lien entre le couple moteur et le couple en sortie de réducteur
1
planétaire s’écrit cr (t) = cm (t).
λ

II.4 Transmission de l’effort à partir de la liaison hélicoïdale

On isole le pignon 2 qui est soumis au couple
( − transmis par le pignon 1 (à travers le contact en I), noté
→ )
−−→ 0
C1→2 = c12 (t)−
→
y0 . On note alors {T1→2 } = −−→ le torseur des actions mécaniques de 1 sur 2. Les mou-
C1→2 J
vements sont suffisamment lents pour qu’on puisse traiter le problème en statique.

Question 7 – Quelles sont les autres actions mécaniques qui s’exercent sur le solide 2 ? Décrire les autres
actions mécaniques par des torseurs.

Correction On isole le pignon 2, il est soumis à :

( −
→ )
0
• l’action du pignon 1 : {T1→2 } = −−→ ;
C1→2 J
X32 −
→x0 + Y32 −
→
y0 + Z32 −→
 
z0
• l’action transmise par 3 à travers la pivot {T3→2 } = ;
L32 −
→
x0 + N32 −
→
z0 J
X02 −
→
x0 + Y02 − →
y0 + Z02 −
→
 
z0
• l’action transmise par 0 à travers la liaison hélicoïdale : {T0→2 } = .
L − →
x +M −
02 0 02 0
→
z +N − →
z
02 0 J
p
Avec M02 = − V Y02 d’après la définition de la liaison hélicoïdale.
2π

DM PSI* 15/34
 Méthodes Numériques

Question 8 – Appliquer le Théorème du Moment Statique en J projeté selon l’axe −

→
y0 et en déduire l’expression
de la force transmise par le solide 0 sur le solide 2.
−
→
Correction On applique le théorème du moment statique en J projeté selon y0 pour profiter du zéro de la
liaison pivot. Tous les torseurs sont déjà au point J donc :

M02 + c12 (t) = 0

2π
Avec la définition de la liaison hélicoïdale : Y02 = c (t) . C’est la force qui permet de mettre en mou-
pV 12
vement de translation le solide 2 par rapport à R0 .

On admet que la force transmise par le solide 0 sur le solide 2 est en fait égale à l’effort de rappel du ressort,
on montrera en deuxième année que l’équation du mouvement se met sous la forme :

dωm
Ieq (t) = cm (t) − Kres Ktrans y(t)
dt
où Ieq l’inertie équivalente ramenée sur l’arbre moteur.

Question 9 – Exprimer l’équation différentielle du mouvement liant le déplacement y(t) au couple cm (t) en
fonction de Ieq , Ktrans , Kres . En déduire la valeur numérique du facteur d’amortissement et conclure quant à
l’amortissement de la réponse indicielle.

Correction On a la relation ẏ(t) = Ktrans ωm (t), on réécrit alors l’équation du mouvement :

Ieq
ÿ(t) + Kres Ktrans y(t) = cm (t)
Ktrans
s
Kres
C’est l’équation d’un oscillateur harmonique forcé non amorti de pulsation propre ω0 = Ktrans .
Ieq

Question 10 – Justifier alors la nécessité de mettre en place une boucle d’asservissement si on souhaite une
commande en position du système.

Correction C’est un système oscillant donc à la limite de l’instabilité (pôles à partie réelle nulle), d’où la
nécessité de mettre en place une boucle d’asservissement pour piloter convenablement le système.

DM PSI* 16/34
 Méthodes Numériques

III Étude de l’effort d’assistance nécessaire au soutien lombaire

Objectif Proposer un modèle de connaissance de l’asservissement en force, le valider par comparaison avec
une mesure sur un banc d’essai et vérifier les performances de l’actionneur linéaire sur un banc d’essai. Ce
modèle permettra de valider une commande pour le cas spécifique étudié.

III.1 Mise en place d’un modèle de connaissance

L’actionneur linéaire placé sur le banc d’essai et sa commande peuvent être modélisés par le schéma-blocs
représenté Figure 5. Le système est commandé en effort mais on observe une deuxième boucle d’asservisse-
ment en interne : la vitesse de rotation du moteur à courant continu est aussi asservi.

Notations et hypothèses :
• la transformée de Laplace de la fonction a(t) est notée A( p) dans le cas général ;
• les conditions de Heaviside sont supposées vérifiées ;
• Fc ( p) représente la consigne en force de l’asservissement de force, dans le domaine de Laplace ;
• F ( p) représente la force développée par l’actionneur linéaire, dans le domaine de Laplace ;
• K1 est le gain du modulateur d’énergie, il est connu.
Les équations modélisant le comportement du moteur électrique (moteur à courant continu) muni d’une
boucle d’asservissement de l’intensité du courant im (t), sont :
• en supposant le temps de réponse de la boucle de courant négligeable,

u I (t) = Rim (t)

• par application des théorèmes généraux de la dynamique appliqués à l’ensemble des solides en mou-
vement,

dωm (t)
Ieq = c m ( t ) − cr ( t ) avec cr (t) = Ty(t)
dt
• loi de couplage électromécanique,

cm (t) = k c im (t)
avec :
• u I (t), la consigne en tension de la boucle d’asservissement de l’intensité du courant im (t) (en V) ;
• im (t), l’intensité du courant d’induit absorbé par le moteur à courant continu (en A) ;
• R, la résistance d’induit du moteur (en Ω) ;
• k c , la constante de couple (en N.m.A−1 ) ;
• Ieq , le moment d’inertie équivalent des solides en mouvement par rapport au référentiel lié au bâti
supposé galiléen, reportée sur l’arbre moteur (en kg.m2 ).
On rappelle qu’on note Ktrans tel que ẏ(t) = Ktrans ωm (t).

Question 11 – Après avoir transformé les équations précédentes dans le domaine de Laplace, exprimer les
gains K3 et K5 en fonction de k c et T.

Correction Dans le domaine de Laplace et à conditions de Heaviside :

• U I ( p) = RIm ( p) ;
• Ieq pΩm ( p) = Cm ( p) − Cr ( p) ;
• Cr ( p) = T · Y ( p) ;
• Cm ( p) = k c Im ( p).
Ainsi on identifie dans le schéma-blocs : K3 = k c et K5 = T .

Question 12 – Exprimer la fonction de transfert H6 ( p) en fonction de Ktrans .

DM PSI* 17/34
 Méthodes Numériques

Ktrans
Correction On a ẏ(t) = Ktrans ωm (t) soit dans le domaine de Laplace Ωm ( p) = Y ( p). On identifie
p
Ktrans
H6 ( p) = .
p

Question 13 – Rappeler sous quelle hypothèse on a le droit d’étudier l’erreur en régime permanent du
système décrit par le modèle de connaissance de la Figure 5.

Correction On peut évaluer l’erreur en régime permanent du système en utilisant le Théorème de la Valeur
Finale (TVF) si et seulement si le système est stable.

Question 14 – En supposant l’hypothèse précédente valide, établir l’expression de Kadapt en fonction des
constantes de la Figure 5 de façon à ce que l’écart ε f (t) en régime permanent soit nul si l’erreur en régime
permanent est nulle.

Correction On ramène la grande boucle de retour (celle avec le gain Kcapt ) sur la sortie F ( p), le gain sur la
Kcapt
boucle de retour devient alors .
Kres
Kcapt
Ainsi pour valider la condition demandée il faut Kadapt = .
Kres

DM PSI* 18/34
 Méthodes Numériques

Figure 5 – Schéma-blocs de l’asservissement de force développée par un actionneur linéaire placé sur le banc
d’essai

DM PSI* 19/34
 Méthodes Numériques

III.2 Réglage de la boucle d’asservissement de la vitesse angulaire du moteur

Le schéma-blocs décrivant la structure de l’asservissement de la vitesse angulaire du moteur est fourni sur la
Figure 6. Cet asservissement doit respecter le cahier des charges fourni dans le Tableau 2.

Figure 6 – Schéma-blocs de la boucle d’asservissement de la vitesse angulaire du moteur électrique

Critère concepteur Niveau

Marge de phase M φ ≥ 80◦
Erreur en régime permanent pour une consigne en échelon constante Nulle
Erreur en régime permanent pour une perturbation en échelon constante Nulle
Pulsation de coupure à 0 dB ω0dB ω0dB = 10 rad/s

Tableau 2 – Critères concepteur pour la boucle d’asservissement de la vitesse angulaire

Le choix d’un correcteur proportionnel intégral est fait par le concepteur. La fonction de transfert du correc-
teur de la boucle d’asservissement en vitesse angulaire est noté Cv ( p), tel que
1 + τi p
Cv ( p) = Ki
τi p
Question 15 – Démontrer qu’il n’y a pas besoin d’utiliser de correcteur pour valider le critère d’erreur en
régime permanent pour une consigne en échelon constante. Quelle opération mathématique (ou quel bloc de
la Figure 6) permet de vérifier à coup sûr ce critère ?

Correction Version 1

En l’absence de correction (Cv ( p) = 1), la fonction de transfert en boucle fermée (sans perturbation puis-
qu’on veut l’erreur et pas la sensibilité à la perturbation) s’écrit :
K1 K3
Ωm ( p) R
HBF ( p) = =
Ωc ( p) K K3
Ieq p + 1
R
C’est un système du premier ordre dont les coefficients au dénominateur sont de même signe, le système
Ω0
est stable. Par TVF et avec Ωc ( p) = on a ε = lim p (Ωc ( p) − Ωm ( p)) = lim Ω0 (1 − HBF ( p)) = 0.
p p →0 p →0

Le critère d’erreur pour une consigne en échelon est validé sans correcteur. En pratique cela est possible
1
grâce à l’intégration (bloc ) dans la boucle ouverte.
p
Version 2

On a une intégration dans la boucle ouverte (BO de classe 1) donc le système est précis pour une consigne
en échelon en absence de perturbation. Ce résultat est vrai sans correcteur.

Question 16 – Démontrer que quelle que soit la valeur de Ki et τi , le correcteur Cv ( p) permet de valider
le critère d’erreur en régime permanent pour une perturbation en échelon. Un correcteur proportionnel
suffirait-il ?

DM PSI* 20/34
 Méthodes Numériques

Correction Cette fois-ci on annule la consigne (Ωc ( p) = 0) et on calcule la fonction de transfert en régula-
C0
tion avec Cr ( p) = :
p
Ωm ( p) 1
Hreg ( p) = =−
Cr ( p) K1 K3 1 + τi p
Ieq p + Ki
R τi p
On calcule la sensibilité à la perturbation s = lim pHreg ( p)Cr ( p) et on voit qu’elle est nulle. Ceci est dû à
p →0
l’intégration en amont de la perturbation en échelon apportée par le correcteur PI.

Sans cette intégration (donc avec un simple correcteur P), on serait sensible à la perturbation en échelon.

Remarque : vous remarquerez que je reprécise toujours la nature des entrées (en échelon ici), cela fait
partie de l’argumentaire pour ce genre de réponse car si on avait une entrée en rampe on n’aurait pas le
même résultat...

Ωm ( p)
On note HBOv ( p) = la fonction de transfert en boucle ouverte de l’asservissement de vitesse angulaire
ε v ( p)
du moteur.

Question 17 – Déterminer l’expression littérale de la phase de HBOv (jω ). En déduire la valeur numérique de
τi respectant les critères concepteur de la boucle de vitesse. On montrera que τi ≥ 0, 57s.

Ki K1 K3 1 + τi p
Correction Dans le domaine de Laplace : HBOv ( p) = · .
Ieq R τi p2
On calcule la phase en degrés dans le domaine fréquentiel : φ(ω ) = Arg ( HBOv (jω )) = arctan(τi ω ) − 180.

On veut une marge de phase M φ = 180 + φ(ω0dB ) ≥ 80◦ , autrement dit arctan(τi ω0dB ) ≥ 80◦ . Par crois-
tan(80)
sance de la fonction arctan on trouve τi ≥ = 0, 57 s .
ω0dB

Le diagramme de Bode de la boucle ouverte HBOv ( p) avec Ki = 1 et τi = 0, 57s déterminé à la question

précédente, est donné sur la Figure 7.

Question 18 – Déterminer la valeur numérique de Ki afin que la boucle d’asservissement de vitesse respecte
les critères concepteur du Tableau 2.

DM PSI* 21/34
 Méthodes Numériques

Figure 7 – Diagramme de Bode de HBOv ( p)

Correction On veut une pulsation de coupure ω0dB = 10 rad/s, autrement dit d’après le diagramme de
gain il faut diminuer le gain de 80dB environ (voir ci-dessous). Alors :

20 log(Ki ) = −80 ⇔ Ki = 10−4 V.s/rad

Tous les critères du tableau 4 sont respectés (marge de phase très légèrement inférieure à celle cherchée
toutefois), le critère de précision et d’insensibilité à la perturbation l’étant forcément par la nature du cor-
recteur PI.

Remarque : le choix de l’unité de Ki est motivé par la dimension (unitaire je pense) du gain K1 .

DM PSI* 22/34
 Méthodes Numériques

III.3 Simplification du modèle de connaissance

Il est possible de mettre le schéma-blocs de la Figure 5 sous la forme du schéma-blocs de la Figure 8, afin de
faciliter la prévision des performances simulées.

Figure 8 – Schéma-blocs de l’asservissement de la force développée par un actionneur linéaire

Question 19 – Déterminer les fonctions de transfert H8 ( p) et H9 ( p) en fonction de K5 , Ieq et H6 ( p). Ne pas

remplacer K5 et H6 ( p) par les expressions trouvées précédemment.

Correction On décale la boucle du capteur angulaire de vitesse d’un cran vers la droite, on trouve alors
1
immédiatement H9 ( p) = .
H6 ( p)
On peut ensuite appliquer la formule de Black à la boucle qui possède la chaîne de retour K5 , on trouve
H6 ( p)
Y ( p) Ieq p H6 ( p)
= , par conséquent H8 ( p) = .
QCm ( p) H6 ( p) H6 ( p)
1+ K5 1+ K5
Ieq p Ieq p

Pour faciliter l’analyse des performances simulées, le schéma-blocs de la Figure 8 est adapté afin de disposer
d’un schéma-blocs à retour unitaire, tel que décrit sur la Figure 9.

Figure 9 – Schéma-blocs de l’asservissement de la force développée par un actionneur linéaire à retour

unitaire

Question 20 – Déterminer l’expression du gain K10 en fonction de Kcapt et de Kres .

Kcapt
Correction On a déjà montré en question 18 que Kadapt = , en factorisant dans la Figure 8 on trouve
Kres
Kcapt
le K10 = de la Figure 9.
Kres
Remarque : l’écart ε f ( p) de la Figure 9 ne peut pas être le même que celui de la Figure 8. En effet dans
Kcapt
la Figure 8 on a ε f ( p) = ( Fc ( p) − F ( p)) et dans la Figure 9 ε f ( p) = Fc ( p) − F ( p). Le ε f ( p) de la
Kres
Figure 9 aurait dû se trouver derrière le gain K10 .

Question 21 – Déterminer la fonction de transfert G ( p) en fonction de H2 ( p), Ieq , H8 ( p), H9 ( p) et Kres . Ne

pas remplacer H2 ( p), H8 ( p) et H9 ( p) par les expressions trouvées précédemment.

DM PSI* 23/34
 Méthodes Numériques

Kres H2 ( p) H8 ( p)
Correction Par formule de Black sur la boucle interne on trouve G ( p) = .
Ieq p + H2 ( p) H8 ( p) H9 ( p)

Pour la suite, on donne la fonction de transfert G ( p), obtenue avec les valeurs de réglage correctes du correc-
teur Cv ( p) :

1 + 0, 57p 1, 2 · 10−5
G ( p) =
p 2 · 10 + 9, 7 · 10−5 p + 5, 3 · 10−6 p2
− 4

III.4 Analyse des performances de l’asservissement en force développée par un action-

neur linéaire
Il est proposé dans cette section d’analyser les performances simulées de l’asservissement en force dont un
extrait du cahier des charges est présenté dans le Tableau 3.

Id Exigence Critère Niveau

Marge de phase Mϕ ≥ 60◦
Id1.1 Stabilité Marge de gain MG > 20 dB
Dépassement maximal < 2, 5%
Id1.2 Précision Erreur en régime permanent pour une entrée en échelon < 1%
Temps de réponse à 5% pour une consigne en échelon de force de 40 N tr5% < 1 s
Id1.3 Rapidité
Vitesse maximale de montée de la force de traction 100 N.s−1

Tableau 3 – Extrait du cahier des charges fonctionnel de l’actionneur linéaire sur le banc d’essai

F ( p)
On note HBO f ( p) = la fonction de transfert en boucle ouverte de l’asservissement en force développé
ε f ( p)
par un actionneur linéaire. Dans un premier temps, le choix d’un correcteur proportionnel C ( p) = Kcor est
F ( p)
réalisé. Le diagramme de Bode de la fonction de transfert HBO f ( p) = = Kcor K10 G ( p) avec Kcor = 1, est
ε f ( p)
donné sur la Figure 10.

Figure 10 – Diagramme de Bode de HBO f ( p)

DM PSI* 24/34
 Méthodes Numériques

Question 22 – Déterminer la valeur numérique limite de Kcor afin que la boucle d’asservissement de force
respecte les critères de marge de phase et de gain du Tableau 3.

Correction La marge de gain est infinie ici car la phase n’atteint jamais −180◦ , regardons la marge de
phase. Elle doit être de 60◦ minimum, autrement dit il faut que la pulsation de coupure se situe à la phase
−120◦ au maximum. Il faut alors augmenter le gain de 75dB minimum pour satisfaire le critère de marge
de phase. Donc :

75
20 log(Kcorr ) ≥ +75 ⇔ Kcorr ≥ 10 20 ≈ 5624rad/(s.V)

Les courbes sur la Figure 11 représentent les réponses temporelles du modèle de connaissance de la Figure 5,
avec les correcteurs Cv ( p) et C ( p) correctement réglés, et de l’expérimentation sur le banc d’essai pour une
consigne en échelon de force de 40 N.

Figure 11 – Réponses temporelles du modèle et expérimentale, pour une consigne en échelon de force de
40N

Question 23 – Quels critères du tableau des exigences (Tableau 3) ne sont pas pris en compte dans le modèle
de connaissance ? D’après la courbe expérimentale, ces critères sont-ils respectés par le système réel ?

Correction Sur la Figure 11 la simulation s’arrête quand le dépassement vaut 2,5%, soit 41N. On ne peut
pas évaluer le temps de réponse du système puisque le régime stationnaire pour une consigne en échelon
n’est pas atteint.

Sur le relevé expérimental ce temps de réponse est clairement inférieur à la valeur 1s demandée par le
cahier des charges et il n’y a pas de dépassement (saturation à 40N). Le cahier des charges concernant le
temps de réponse à 5% est validé expérimentalement.
36 − 4
La vitesse de montée expérimentale vaut environ = 80N/s < 100N/s donc le critère de vitesse
0, 5 − 0, 1
de montée est aussi respecté. Toutefois la vitesse de montée maximale relevée sur la courbe de simulation
41 − 7
vaut environ = 212N/s > 100N/s, le critère de vitesse en montée n’est pas respecté lors de la
0, 28 − 0, 12

DM PSI* 25/34
 Méthodes Numériques

simulation.

IV Amélioration du modèle. Mise en place d’une limitation en vitesse

angulaire
Pour améliorer le modèle de connaissance et le valider, la comparaison entre la réponse simulée issue du
modèle de connaissance amélioré et la réponse expérimentale sera traitée par résolution numérique informa-
tique. Le langage de programmation utilisé est Python.

Figure 12 – Schéma-blocs de l’asservissement de force développée par l’actionneur linéaire avec limitation
de la vitesse angulaire

Notations et hypothèses :
• f (t) est la grandeur de sortie de l’asservissement en force, de variable de Laplace F ( p) ;
• f c (t) est la grandeur de consigne de l’asservissement en force, de variable de Laplace Fc ( p) ;
• ωc (t) est la commande de vitesse angulaire du moteur, de variable de Laplace Ωc ( p) ;
• la dérivée première temporelle d’une fonction h(t) est notée ḣ(t) et sa dérivée seconde ḧ(t) ;
• les conditions de Heaviside sont supposées vérifiées, soient f (t = 0) = 0 et f˙(t = 0) = 0.

IV.1 Cadre général de la résolution numérique d’un problème de Cauchy

Pour déterminer numériquement la solution E(t) d’une équation différentielle, il faut préalablement la mettre
sous la forme d’un problème de Cauchy :

Ė(t) = FCauchy ( E(t), t) avec E(t = 0) = E0 (1)

• E(t) est appelé vecteur d’états à l’instant t ;
• Ė(t) représente le vecteur composé des dérivées premières des états par rapport au temps t ;
• FCauchy est la fonction de Cauchy ;
• E0 est le vecteur d’états des conditions initiales.
Le problème étant décrit sous la forme de Cauchy, la résolution numérique peut être menée par un schéma
numérique du type Euler explicite ou en utilisant la fonction odeint de la bibliothèque scipy. C’est un
schéma d’Euler explicite qui va être mis en place dans ce qui suit.

DM PSI* 26/34
 Méthodes Numériques

IV.2 Étapes de la résolution du problème

Dans la Figure 12, on remarque une zone de non linéarité. On propose de traiter celle-ci à part. De plus
compte-tenu du correcteur Cv ( p) de la boucle d’asservissement de la vitesse angulaire du moteur électrique
réglé aux questions précédentes, il est possible de simplifier la fonction de transfert G ( p) qui s’écrit mainte-
nant :
1, 18
G ( p) =
p(18 + p)
Ainsi on décompose la résolution du problème de la manière suivante :
1. écriture de la fonction de Cauchy associée à la fonction de transfert G ( p), notée G(E,wc,t) ;
2. écriture de la fonction de Cauchy pour la boucle fermée de l’asservissement en effort sans non linéarité
(pour s’échauffer), notée FTBFsnl(E,t) ;
3. écriture de la fonction de Cauchy pour la boucle fermée de l’asservissement en effort avec non linéarité
(pour se muscler), notée FTBFanl(E,t) ;
4. écriture d’une fonction de tracé de la réponse temporelle obtenue pour le modèle de la Figure 12 avec
non linéarité, notée Trace(E0,t_fin,N_pt).
Les notations sont les suivantes :
• E est une liste de 2 nombres représentant le vecteur d’états ;
• E0 est une liste de 2 nombres représentant le vecteur d’états au temps t = 0 ;
• wc est un nombre représentant la consigne de vitesse angulaire du moteur ;
• t est un nombre représentant le temps ;
• t_fin est la durée de l’essai ;
• N_pt est le nombre de points souhaité sur le graphe.

IV.2.1 Fonction de Cauchy associée à la fonction de transfert G ( p)

f˙(t)
   
f (t)
Le vecteur d’états associé à l’étude de G ( p) retenu est E(t) = ˙ , ainsi Ė(t) = ¨ = G ( E ( t ), ω c ( t ), t ).
f (t) f (t)
Question 24 – Déterminer la fonction G ( E(t), ωc (t), t), associée à G ( p).

Correction
f˙(t) f˙(t)
     
E[1](t)
On a Ė(t) = = = avec E[1](t) = f˙(t) la
f¨(t) −18 f˙(t) + 1, 18ωc (t)
−18E[1](t) + 1, 18ωc (t)
 
E[1](t)
deuxième coordonnée du vecteur E(t). Alors G ( E(t), ωc (t), t) = .
−18E[1](t) + 1, 18ωc (t)

Question 25 – Écrire en langage Python la fonction G(E,wc,t) qui implémente la fonction G ( E(t), ωc (t), t).
Cette fonction renvoie une liste de 2 nombres à virgule flottante.

Correction

1 def G(E,wc ,t):

2 G0 = E[1] # Première coordonnée
3 G1 = −18 * E[1] + 1.18 * wc # Deuxième coordonnée
4 r e t u r n [G0 ,G1]

DM PSI* 27/34
 Méthodes Numériques

IV.2.2 Fonction de Cauchy pour la boucle fermée de l’asservissement en effort sans non linéarité
Le correcteur de la boucle d’asservissement de la force développée par un actionneur linéaire est un correcteur
proportionnel tel que C ( p) = Kcor . Ainsi, en l’absence de limitation, la grandeur ωc (t) est telle que :

ωc (t) = Kcor K10 ε f (t) = Kcor K10 ( f c (t) − f (t))

Il reste à exprimer la fonction de Cauchy pour la boucle fermée de l’asservissement de force :

f˙(t)
 
Ė(t) = ¨ = FTBFsnl ( E(t), t)
f (t)
On considère dès à présent que la consigne de force f c (t) est un échelon d’amplitude Fcons , prise égale à 40 N.

Question 26 – Recopier et compléter le programme Python ci-dessous implémentant la fonction FBFsnl qui
renvoie une liste de deux nombres.
1 def FTBFsnl (E,t):
2 Fcons = 40 # échelon de 40N
3 K10 = 0.0277 # gain d’adaptation
4 Kcor = 5400 # gain du correcteur proportionnel
5 wc = # à compléter
6 # compléter en terminant la fonction

Correction

1 def FTBFsnl (E,t):

2 Fcons = 40 # échelon de 40N
3 K10 = 0.0277 # gain d’adaptation
4 Kcor = 5400 # gain du correcteur proportionnel
5 wc = K10 * Kcor * ( Fcons − E[0]) # Vitesse angulaire de consigne
6 r e t u r n G(E,wc ,t)

IV.2.3 Fonction de Cauchy pour la boucle fermée de l’asservissement en effort avec non linéarité
Question 27 – Recopier et compléter le programme Python ci-dessous implémentant la fonction FBFanl qui
renvoie une liste de deux nombres et prend en compte la limitation de la consigne de vitesse angulaire du
moteur.
1 def FTBFanl (E,t):
2 Fcons = 40 # échelon de 40N
3 K10 = 0.0277 # gain d’adaptation
4 Kcor = 5400 # gain du correcteur proportionnel
5 wcmax = 1250 # limitation en vitesse angulaire dans l’intervalle [-wcmax,wcmax]
6 # compléter en terminant la fonction

Correction

1 def FTBFanl (E,t):

2 Fcons = 40 # échelon de 40N
3 K10 = 0.0277 # gain d’adaptation
4 Kcor = 5400 # gain du correcteur proportionnel
5 wcmax = 1250 # limitation en vitesse angulaire dans l’intervalle [-wcmax,wcmax]
6 wc = K10 * Kcor * ( Fcons − E[0]) # Vitesse angulaire de consigne
7 # Prise en compte de la saturation
8 i f wc > wcmax :
9 wc = wcmax
10 i f wc < −wcmax :

DM PSI* 28/34
 Méthodes Numériques

11 wc = −wcmax
12 r e t u r n G(E,wc ,t)

IV.2.4 Tracé de la simulation et conclusion vis-à-vis du cahier des charges

Question 28 – Quel autre nom donne-t-on à la limitation de vitesse angulaire : saturation ou seuil ?

Correction C’est une saturation en vitesse.

Les scripts précédemment écrits ainsi qu’une fonction de tracé ont permis d’obtenir la courbe de la Figure 13.

Question 29 – Écrire la fonction Trace(E0,t_fin,N_pt) qui prend en arguments des conditions initiales
E0, une durée t_fin (départ à t = 0) et un nombre de points N_pt et qui renvoie la courbe de simulation.
Celle-ci doit être la même que celle tracée en Figure 13 pour vérifier votre résultat.

Correction

1 # Tracé de la courbe de simulation

2 def Trace (E0 ,t_fin ,N_pt):
3 # Construction de la liste des temps
4 delta_t = t_fin /( N_pt −1) # On a N_pt-1 intervalles de temps
5 liste_t = [k * delta_t f o r k in range (N_pt)]
6
7 # Initialisation de la liste des vecteurs d’états E aux différents instants
8 liste_E = [E0]
9
10 # Résolution avec un schéma d’Euler explicite
11 f o r i in range (1, N_pt):
12 E_new0 = liste_E [i −1][0] + delta_t * FTBFanl ( liste_E [i−1], liste_t [i
−1]) [0]
13 E_new1 = liste_E [i −1][1] + delta_t * FTBFanl ( liste_E [i−1], liste_t [i
−1]) [1]
14 liste_E . append ([ E_new0 , E_new1 ])
15
16 # Liste des valeurs de force calculées à chaque instant
17 liste_f = [ liste_E [k][0] f o r k in range (N_pt)]
18
19 # Tracé de la courbe
20 fig = plt. figure ()
21 [Link](liste_t , liste_f )
22 [Link](’both ’)
23 plt. xlabel (’temps (s)’)
24 plt. ylabel (’force (N)’)
25 [Link] ()
26
27 # Tracé de la réponse à conditions initiales nulles
28 Trace ([0 ,0] ,1 ,200)

Question 30 – Conclure sur la capacité du correcteur proportionnel (de gain 5400) avec limitation de la vitesse
angulaire à respecter le cahier des charges, en analysant l’écart entre les performances simulées (Figure 13)
et les performances attendues. Se limiter aux critères pertinents du Tableau 3.

DM PSI* 29/34
 Méthodes Numériques

Correction Le système est stable, ensuite :

41 − 40
• on a une valeur maximale de 41N environ, soit un dépassement de = 2, 5% (dépassement au
40
maximum autorisé) ;
• le système est précis pour une entrée en échelon (critère Id1.2 validé) ;
• le temps de réponse à 5% est inférieur à 0,6s (< 1s), le critère de rapidité sur le temps de réponse est
validé ;
36 − 4
• la pente maximale relevée sur la courbe vaut = 80N/s < 100N/s, le critère de vitesse de
0, 5 − 0, 1
montée est respecté.
Ainsi tous les items du cahier des charges sont validés.

Figure 13 – Réponse temporelle du modèle de connaissance pour une consigne de 40 N avec correction
proportionnelle et limitation de la vitesse angulaire

V Synthèse et ouverture de l’étude

Objectif Valider le modèle de connaissance, valider la commande optimisée et envisager un prolongement
à l’étude.

Remarque Attention, les questions qui vont suivre demandent des réponses qualitatives, il faut être très
vigilant à répondre en chiffrant vos propos.

V.1 Validation du modèle de connaissance et de la commande optimisée

On appelle commande optimisée la commande avec le correcteur de la boucle d’asservissement de force du
type proportionnel avec limitation de la vitesse angulaire. Sur le banc d’essai avec l’actionneur linéaire, on
implante ce correcteur, et on procède à quatre essais pour des consignes en échelon de force de 10, 20, 30 et
40 N. Les mesures correspondantes sont tracées en traits discontinus sur la Figure 14. Ce même correcteur
est mis en place dans le modèle de connaissance de la Figure 5, on obtient les réponses temporelles de
simulation tracées en traits continus sur la Figure 14.

DM PSI* 30/34
 Méthodes Numériques

Figure 14 – Résultats des simulations et des expérimentations pour une entrée de consigne de force en
échelon d’amplitude 10, 20, 30 et 40 N

La démarche de l’ingénieur abordée dans le programme de sciences industrielles de l’ingénieur de la filière

MP et PSI s’appuie sur les écarts entre trois performances (Figure 15).

Figure 15 – Synoptique de la démarche de l’ingénieur, telle que présentée dans le programme

En se référant à la Figure 14, l’analyse de l’écart entre les performances simulées et les performances mesu-
rées valide le modèle de connaissance de l’asservissement de force.

Question 31 – Choisir un des écarts L-C, S-L ou S-C permettant de valider la commande optimisée. Effectuer
l’analyse de cet écart. Il est attendu une argumentation rigoureuse s’appuyant sur les données et les références
du texte. Les numéros de figure, de tableau, ou d’exigence sont, par exemple, des références utilisables.

Correction Dans ce sujet, un banc d’essais a été monté pour réaliser des expériences et une modélisation
(figure 11) basée sur ce banc d’essais a été proposée afin de prévoir les performances du système réel. Ainsi
c’est l’écart S-L qui doit être analysé.

Sur la Figure 14 on observe :

• Consigne de 10N : il n’y a pas de dépassement ni en simulation ni en expérimentation. L’erreur
statique est nulle dans les 2 cas. La vitesse de montée est plus lente lors de l’expérimentation et le

DM PSI* 31/34
 Méthodes Numériques

temps de réponse à 5% est plus grand (0,25s contre moins de 0,2s pour la simulation). Les frottements
secs ne sont donc pas négligeables pour cette consigne, or ils n’ont pas été pris en compte dans la
modélisation de la figure 11. Les exigences de dépassement, rapidité et précision du tableau 5 sont
respectées.
• Consigne de 20N : un très léger dépassement est visible en simulation mais inférieur à 2,5% (maxi-
mum autorisé par l’exigence de stabilité du tableau 5), il n’y a pas de dépassement lors de l’ex-
périence. Les observations sont similaires à celles faites pour le cas de la consigne de 10N mais le
système réel met toujours un peu plus de temps à atteindre sa vitesse de montée maximale (les
courbes ne se superposent bien qu’au delà de 1,15s sur la Figure 14). Le frottement sec joue encore
pour cette valeur de consigne.
• Consignes de 30N et 40N : cette fois les courbes expérimentales et de simulations se superposent bien
au bruit de mesure près. Le frottement sec a un effet négligeable et toutes les exigences du tableau 5
sont respectées.
Globalement la modélisation est fidèle au système réel. Afin d’améliorer le modèle on pourrait envisager
de prendre en compte les effets des frottements présents dans le banc d’essais.

V.2 Étude du système perturbé

Le schéma-blocs de la Figure 18 introduit une perturbation Ypert ( p), représentant une perturbation dans le
domaine temporel ypert (t).

Question 32 – Quelle est l’unité de ypert (t) ? Que peut représenter cette perturbation dans le contexte de
l’exosquelette lombaire ? En se référant à la problématique du sujet, quel est l’intérêt d’introduire cette per-
turbation dans le modèle ?

Correction ypert (t) est homogène à une longueur en mètres. Dans le cadre de l’exosquelette lombaire, qui
s’assimile à 2 ceintures liées par des actionneurs verticaux, on peut très bien imaginer des glissements entre
la ceinture et l’utilisateur ce qui peut générer des perturbations de déplacement. Il est donc nécessaire de
corriger les éventuels défauts de positionnement des ceintures afin de maintenir une pression lombaire au
niveau prévu par le commanditaire.

On effectue une simulation avec une perturbation normalisée représentée sur la Figure 16.

Figure 16 – Évolution temporelle de la perturbation ypert (t)

La Figure 17 montre l’évolution de la force développée par un actionneur linéaire soumis à la perturbation
décrite sur la Figure 16. C’est le résultat d’une simulation obtenue à partir du modèle de connaissance avec
le correcteur proportionnel et limitation sur la vitesse angulaire.

DM PSI* 32/34
 Méthodes Numériques

Figure 17 – Évolution de la force développée par l’actionneur linéaire soumis à une perturbation. Consigne
en échelon de 40 N avec perturbation à t = 0, 8 s

Question 33 – Analyser la courbe de simulation de la Figure 17 et conclure au regard de la problématique

du sujet.

Correction La perturbation est une fonction rampe qui démarre à partir de t = 0, 8s et devient constante
à partir de t = 2s. Si on note u(t) la fonction de Heaviside on a : ypert (t) = Y0 (t − 0, 8)u(t − 0, 8) − Y0 (t −
2) u ( t − 2).
1 1
Dans le domaine de Laplace la perturbation est en sur la partie rampe et en à partir de t = 0, 2s. Or
p2 p
dans le schéma-blocs de la Figure 18 on a une double intégration en amont de la perturbation grâce aux
1
blocs et au bloc H6 ( p) identifié en question 16. Le système est alors insensible à une perturbation en
Ieq p
rampe et à une perturbation constante.

Cela explique les 2 "bosses" qui s’écrasent sur la Figure 17.

Question 34 – Le banc d’essai équipé d’un actionneur linéaire, dans la configuration étudiée dans ce sujet,
permet-il d’analyser l’écart S-L en présence d’une perturbation ypert (t) ?

Correction On a déjà étudié l’écart S-L en question 34 puisque le modèle se basait sur le système du
laboratoire. Toutefois le système du laboratoire est à l’horizontale et on a supposé en question 35 que la
perturbation venait de défauts de positionnement de l’exosquelette sur l’utilisateur (donc avec les action-
neurs à la verticale).

Enfin la proposition de perturbation de la Figure 16 n’est peut-être pas réaliste. Il n’est pas possible
d’évaluer un écart S-L quant à un défaut de positionnement.

DM PSI* 33/34
 Méthodes Numériques

Figure 18 – Schéma-blocs de l’asservissement de force développée par un actionneur linéaire placé sur le
banc d’essai avec perturbation

DM PSI* 34/34