6. Contrôle de trajectoires#
Marc BUFFAT, dpt mécanique, Université Claude Bernard Lyon 1
Mise à disposition selon les termes de la Licence Creative Commons Attribution - Pas d’Utilisation Commerciale - Partage dans les Mêmes Conditions 2.0 France
%matplotlib inline
import numpy as np
import sympy as sp
import matplotlib.pyplot as plt
from IPython.core.display import HTML
from IPython.display import display,Image
from sympy.physics.vector import init_vprinting
init_vprinting(use_latex='mathjax', pretty_print=False)
plt.rc('font', family='serif', size='14')
6.1. Problème: pont roulant#
6.1.1. objectif#
On veut déplacer, à l’aide d’un système de levage par câble (grue ou pont roulant), une charge d’un endroit à un autre en évitant les oscillations résiduelles à l’arrivée. L’objectif est de proposer une stratégie de commande simple et réaliste qui repose sur une structure de commande hiérarchisée, composée de régulateurs de bas niveau rapides, simples et découplés et d’une commande de haut niveau lente et prenant en compte les couplages. En outre, on mesure la position et la vitesse du chariot ainsi que la position et la vitesse du treuil, mais la position de la charge n’est pas mesurée. C’est pourquoi nous ne considérons que des bouclages qui ne dépendent que des positions et vitesses du chariot et du treuil.
6.1.2. modélisation#
modèle mécanique;
pendule P de masse m de longueur \(l(t)\) accroché à un chariot C de masse M se deplacant horizontalement \(x_c(t)\)
système à 3 ddl \(xc(t), \theta(t), l(t) \)
6.1.3. repère et ddl#
from sympy.physics.mechanics import dynamicsymbols, Point, ReferenceFrame
from sympy.physics.mechanics import Particle, RigidBody, Lagrangian, LagrangesMethod
# parametres
theta, xc, l = dynamicsymbols("theta x_c l")
thetap, xcp, lp = dynamicsymbols("theta x_c l",1)
m,M,g = sp.symbols('m M g')
t = sp.Symbol('t')
# pt et répère
O = Point('O')
RO = ReferenceFrame("R_O")
C = Point('C')
C.set_pos(O,xc*RO.x)
RP = ReferenceFrame("R_P")
RP.orient(RO,"Axis",[theta,RO.z])
P = Point('P')
P.set_pos(C, -l*RP.y)
display(C.pos_from(O))
display(P.pos_from(O))
6.1.4. cinématique#
O.set_vel(RO,0.)
C.set_vel(RO,xcp*RO.x)
P.set_vel(RP, -lp*RP.y)
display(C.vel(RO))
display(P.vel(RP))
P.v1pt_theory(C,RO,RP)
Pc = Particle('P_c',C,M)
Pa = Particle('P_a',P,m)
display(Pa.kinetic_energy(RO))
Pa.potential_energy = m*g*Pa.point.pos_from(C).dot(RO.y)
display(Pa.potential_energy)
6.1.5. Lagrangien#
La = Lagrangian(RO,Pa,Pc)
display(La)
6.1.6. Bilan des forces#
force motrice \(F_c\) sur le chariot en C
force de tension \(T\) dans le cable (ne travaille pas)
couple sur un treuil de rayon \(\rho\)
rayon \(\rho\) (on néglige son inertie \(J_\rho\))
angle \(\phi = (l-l_0)/\rho\)
couple \(C_t\)
travail \(C_t \rho \delta l \)
Fc, Ct, rho, l0 = sp.symbols('F_c C_t rho l_0')
phi_t = (l-l0)/rho
Rt = ReferenceFrame('R_t')
Rt.orient(RO,'Axis',[phi_t,RO.z])
display(La)
LM = LagrangesMethod(La,[theta,xc,l],forcelist=[(C,Fc*RO.x),(Rt,-Ct*m*rho*RO.z)], frame = RO)
eq = sp.simplify(LM.form_lagranges_equations())
display(eq)
6.2. Chariot immobile#
\(\dot{x_c}=0\)
la force de tension ne travaille pas (liaison parfaite)
mouvement \(l(t)\) imposé
LLa = sp.simplify(La.subs({xcp:0}))
display(LLa)
LM = LagrangesMethod(LLa,[theta])
eq=sp.simplify(LM.form_lagranges_equations()[0]/(m*l))
display(eq)
6.2.1. interprétation#
variation moment cinétique suivant Ro.y = moment poids
pz = Pa.angular_momentum(C,RO).subs({xcp:0}).dot(RO.z)
display(pz)
pz.diff(t) + -m*g*(RO.y.cross(P.pos_from(C)).dot(RO.z))
6.2.2. linéarisation#
\(l \approx l_0\)
\(\dot{l} = \omega l_0\)
solution amortie ou amplifiée
l0 , omega = sp.symbols('l_0 omega')
eq0=sp.simplify(eq.subs({lp:omega*l0,l:l0,g:omega**2*l0,sp.sin(theta):theta})/l0)
display(eq0)
theta0 = sp.Symbol('theta_0')
sol = sp.dsolve(eq0,ics={theta.subs(t,0):theta0, theta.diff(t).subs(t,0):0})
display(sol)
6.2.3. simulation cas linéarisé#
Theta=sp.lambdify([t,theta0,omega],sol.rhs,'numpy')
L0 = 10.0
GG = 9.81
OMEGA0 = 0.2*np.sqrt(GG/L0)
THETA0 = np.pi/8
TF = 1.5*np.pi/OMEGA0
tt = np.linspace(0,TF,100)
THETA = Theta(tt,THETA0,OMEGA0)
LL = L0+OMEGA0*L0*tt
plt.subplot(1,2,1)
plt.plot(tt,THETA)
plt.title('angle $\\theta$')
plt.subplot(1,2,2)
plt.plot(tt,LL)
plt.title("allongement l(t)")
Text(0.5, 1.0, 'allongement l(t)')

6.2.4. animation cas linéaire#
from validation.Chariot import Pendule,Chariot
pas = 2
pendule = Pendule(LL[0::pas],THETA[0::pas],tt[pas])
anim = pendule.calculanim(-20,20,-60,1);

plt.rc('animation', html='jshtml')
anim
6.2.5. simulation cas non linéaire#
resolution \(Y = [\theta, \dot\theta, l]\)
equation sur \(l\)
si \(t > T_1\) $\(\frac{dl}{dt} = 0 \)$
si \(t > T_1\) et \(t<2 T_1\) $\(\frac{dl}{dt} = \sqrt{\frac{g}{l}} \)$
sinon $\(\frac{dl}{dt} = 0 \)$
6.2.6. equations de Lagrange à l’ordre 1#
smb = sp.simplify(LM.rhs())
display(smb)
smbF = sp.lambdify([theta,thetap,l,lp,g],smb,'numpy')
6.2.7. simulation numérique#
# parametres
OMEGA0 = np.sqrt(GG/L0)
T1 = 4*2*np.pi/OMEGA0
def F(Y,t):
global L0
Lp = 0.0
if (t>T1) and (t<2*T1): Lp = np.sqrt(GG/Y[2])
FF = smbF(Y[0],Y[1],Y[2],Lp,GG)
dY = np.array([FF[0,0],FF[1,0],Lp])
return dY
from scipy.integrate import odeint
Y0 = [THETA0,0.,L0]
TT = 4*T1
tt = np.linspace(0,TT,200)
sol = odeint(F,Y0,tt)
THETA = sol[:,0]
LL = sol[:,2]
plt.subplot(1,2,1)
plt.plot(tt,THETA)
plt.subplot(1,2,2)
plt.plot(tt,LL)
[<matplotlib.lines.Line2D at 0x7fe896eec700>]

6.2.8. animation#
pas = 2
pendule = Pendule(LL[0::pas],THETA[0::pas],tt[pas])
anim = pendule.calculanim(-10,10,-30,1);

plt.rc('animation', html='jshtml')
anim
6.2.9. Conclusion#
on ne contrôle pas le mouvement de la charge !!
6.3. Chariot mobile, longueur fixe#
on applique une force horizontale \(F_c\) avec un cable de longueur \(l_0\)
adim: \( F_c M = F_c \alpha m \)
\(\alpha = \frac{M}{m}\)
objectif:
déterminer la force \(F_c\) pour avoir un ballant minimum
tests de différents contrôle
contrôle du chariot seule
asservissement de la trajectoire du chariot
couplage de haut niveau
alpha = sp.Symbol('alpha')
LLa = sp.simplify(La.subs({lp:0,l:l0,M:alpha*m}))
display(LLa)
6.3.1. equations de Lagrange#
Fc = sp.symbols('F_c')
LM = LagrangesMethod(LLa,[theta,xc],forcelist=[(C,Fc*alpha*m*RO.x)], frame = RO)
eq = sp.simplify(LM.form_lagranges_equations())
display(eq)
display(LM.mass_matrix)
display(sp.simplify(LM.forcing))
display(sp.simplify(LM.rhs()).doit())
6.3.2. Contrôle simple du chariot#
On contrôle uniquement le chariot \(x_c(t)\)
on impose une trajectoire au chariot tq. $\( x_c(0)=0, \dot{x_c}(0)=0, x_c(t_f) = x_f, \dot{x_c}(t_f)=0\)$
\(x_c(t)\) polynôme de degré 3
d’où la force \(F_c\) sur le chariot
xf,tf = sp.symbols('x_f t_f')
xcc = xf*t**2*(3*tf-2*t)/tf**3
display(xcc)
xcc.diff(t).subs(t,tf)
Fcc = (1+alpha)/alpha*xcc.diff(t,2)
display(Fcc)
# force sur le chariot
FCX = sp.lambdify([t,tf,xf,alpha],Fcc,'numpy')
6.3.3. Asservissement simple#
on applique un asservissement simple pour avoir la trajectoire précédente
asservissement $\( F_a(t) = -K_c \left( (x(t)-x_i(t)) + T_c(\dot x - \dot x_i) \right) \)$
Kc,Tc = sp.symbols('K_c T_c')
Fca = -Kc*((xc-xcc)+Tc*(xcp-xcc.diff(t)))
display(Fca)
# force sur le treuil
FCA = sp.lambdify([t,tf,xf,xc,xcp,Kc,Tc],Fca,'numpy')
6.3.4. Commande de haut niveau lente#
linéarisation des équations
controle du positionnement de la charge
6.3.4.1. mouvement de la charge#
coordonnees \(\xi(t)\), \(\eta(t)\) dans le repere fixe
display(P.pos_from(O).express(RO))
xi, eta = dynamicsymbols('xi eta')
display(sp.Eq(xi,P.pos_from(O).dot(RO.x)))
display(sp.Eq(eta,P.pos_from(O).dot(RO.y)))
display(sp.Eq(sp.tan(theta),(-xi+xc)/eta))
6.3.4.2. parametrisation / trajectoire#
equation d’équilibre de la charge en fonction de la tension T
T = sp.Symbol('T')
display(sp.Eq(m*xi.diff(t,t) , -T*sp.sin(theta)))
display(sp.Eq(m*eta.diff(t,t), T*sp.cos(theta) - m*g))
display(sp.Eq(sp.tan(theta), -xi.diff(t,t)/(eta.diff(t,t)+g)))
display(sp.Eq(xc, xi - eta*xi.diff(t,t)/(eta.diff(t,t)+g)) )
6.3.4.3. paramétrisation#
si on connait \(xi(t)\) et \(eta(t)\) jusqu’à leurs dérivées d’ordre 4, peut donc déterminer le mouvement et calculer la force \(F_c\)
6.3.4.4. linéarisation \(\theta\) petit#
\(\eta = - l_0\)
\(\theta = -\frac{\ddot{\xi}}{g}\)
\(x_c = \xi + \frac{l_0}{g} \ddot{\xi}\)
\(\omega_0 ^2 = \frac{g}{l_0}\)
calcul de la force \(F_c\)
omega0 = sp.Symbol('omega_0')
eq0=sp.simplify(eq[1].subs({sp.cos(theta):1,sp.sin(theta):0})/m)
display(eq0)
eq1=eq0.subs({xc.diff(t,2):xi.diff(t,2)+l0/g*xi.diff(t,4),theta.diff(t,2):-xi.diff(t,4)/g, g: omega0**2*l0})
display(eq1)
Fc1=(alpha*Fc+eq1)/alpha
display(sp.Eq(Fc,Fc1))
6.3.4.5. trajectoire#
trajectoire \(\xi(t)\) tq. $\(\xi(0)=0, \frac{d\xi}{dt}(0)=0, \frac{d^2\xi}{dt^2}(0)=0, \frac{d^3\xi}{dt^3}(0)=0, \frac{d^4\xi}{dt^4}(0)=0\)\( \)\(\xi(t_f)=x_f, \frac{d\xi}{dt}(t_f)=0, \frac{d^2\xi}{dt^2}(t_f)=0, \frac{d^3\xi}{dt^3}(t_f)=0, \frac{d^4\xi}{dt^4}(t_f)=0\)$
\(\xi(t)\) polynome de degré 9 en t
xf, tf = sp.symbols('x_f t_f')
x = xf*(t/tf)**5*(126-420*(t/tf)+540*(t/tf)**2-315*(t/tf)**3+70*(t/tf)**4)
display(sp.Eq(xi,x))
# verification
print("dérivées en tf:",x.diff(t,1).subs({t:tf}),x.diff(t,2).subs({t:tf}),x.diff(t,3).subs({t:tf}),x.diff(t,4).subs({t:tf}))
# contrôle Fc
Fcx=Fc1.subs({xi.diff(t,4):x.diff(t,4),xi.diff(t,2):x.diff(t,2)}).doit().simplify()
display(sp.Eq(Fc,Fcx))
# position chariot xc
xcx= xi + xi.diff(t,2)/omega0**2
xcx= xcx.subs({xi:x}).doit().simplify()
display(sp.Eq(xc,xcx))
# angle theta
thx= -xi.diff(t,2)/(l0*omega0**2)
thx= thx.subs({xi:x}).doit().simplify()
display(sp.Eq(theta,thx))
dérivées en tf: 0 0 0 0
# fonctions numériques
XI = sp.lambdify([t,xf,tf],x,'numpy')
XC = sp.lambdify([t,xf,tf,omega0],xcx,'numpy')
FX = sp.lambdify([t,xf,tf,omega0,alpha],Fcx,'numpy')
THX= sp.lambdify([t,xf,tf,omega0,l0],thx,'numpy')
6.3.5. Parametre de simulation#
attention valeur des parametres KC et TC
optimal si \(Tc^2 = 4/K_c\)
# PARAMETRES NUMERIQUES
# deplacement lent, rapide, + rapide
TF = 8.0
TF = 6.0
TF = 4.0
#
XF = 4.0
ALPHA = 0.5
#ALPHA = 10.0
L0 = 3.0
GG = 9.81
OMEGA0 = np.sqrt(GG/L0)
KC = 100.0
TC = np.sqrt(4.0/KC)
print("Simulation K,TC,TF: ",KC,TC,TF,2*np.pi/OMEGA0)
Nt = 400
tt = np.linspace(0,TF,Nt)
Simulation K,TC,TF: 100.0 0.2 4.0 3.4746094143618937
6.3.6. solution planifiée optimale#
XXI = XI(tt,XF,TF)
XCX = XC(tt,XF,TF,OMEGA0)
FFX = FX(tt,XF,TF,OMEGA0,ALPHA)
TTHX= THX(tt,XF,TF,OMEGA0,L0)
plt.figure(figsize=(10,8))
plt.subplot(3,1,1)
plt.plot(tt,XXI,label="P")
plt.plot(tt,XCX,'--',label="C")
plt.legend()
plt.title("Solution planifiée position P et C")
plt.subplot(3,1,2)
plt.plot(tt,FFX)
plt.ylabel("Force F")
plt.subplot(3,1,3)
plt.plot(tt,TTHX,'--')
plt.ylabel("angle $\\theta$")
plt.xlabel('t');

# tracer du mouvement optimal
plt.figure(figsize=(12,8))
for i in range(0,len(tt),20):
X,Y = [XCX[i],XXI[i]],[0.,-L0]
plt.plot(X,Y,'-o',ms=20)
plt.axis('equal');

6.3.7. Simulation#
display(sp.simplify(LM.rhs()).doit())
smb = sp.simplify(sp.simplify(LM.rhs()).doit().subs({l:l0,g:omega0**2*l0}))
display(smb)
smbF = sp.lambdify([theta,xc,thetap,xcp,Fc,l0,omega0,alpha],smb,'numpy')
def Fs(Y,t):
"""commande simple"""
global L0,OMEGA0,ALPHA,TF,XF
FC = FCX(t,TF,XF,ALPHA)
FF = smbF(Y[0],Y[1],Y[2],Y[3],FC,L0,OMEGA0,ALPHA)
return FF[:,0]
def Fa(Y,t):
"""asservissement simple"""
global L0,OMEGA0,ALPHA,TF,XF,KC,TC
FC = FCA(t,TF,XF,Y[1],Y[3],KC,TC)
FF = smbF(Y[0],Y[1],Y[2],Y[3],FC,L0,OMEGA0,ALPHA)
return FF[:,0]
def F(Y,t):
"""commande lente couplée"""
global L0,OMEGA0,ALPHA
FC = FX(t,XF,TF,OMEGA0,ALPHA)
FF = smbF(Y[0],Y[1],Y[2],Y[3],FC,L0,OMEGA0,ALPHA)
return FF[:,0]
from scipy.integrate import odeint
Y0 = [0.0, 0.0, 0., 0. ]
# perturbation
Y0 = [0.01,0.01, 0.,0.]
tt = np.linspace(0,TF,Nt)
sol = odeint(F,Y0,tt)
sols= odeint(Fs,Y0,tt)
sola= odeint(Fa,Y0,tt)
THETA = sol[:,0]
XXC = sol[:,1]
print("Erreur sur la position finale")
print("commande couplée erreur xc={:.4f} \t theta={:.3f} deg".format(sol[-1,1]-XF ,np.degrees(sol[-1,0])))
print("commande simple erreur xc={:.4f} \t theta={:.3f} deg".format(sols[-1,1]-XF,np.degrees(sols[-1,0])))
print("commande asservie erreur xc={:.4f} \t theta={:.3f} deg".format(sola[-1,1]-XF,np.degrees(sola[-1,0])))
Erreur sur la position finale
commande couplée erreur xc=0.0525 theta=-0.645 deg
commande simple erreur xc=-0.5876 theta=17.988 deg
commande asservie erreur xc=0.0652 theta=13.829 deg
plt.figure(figsize=(12,8))
plt.subplot(2,1,1)
plt.plot(tt,THETA,label="couplée")
plt.plot(tt,sols[:,0],label="simple")
plt.plot(tt,sola[:,0],label="asservie")
plt.plot(tt,TTHX,'--',label="optimale")
plt.ylabel('angle $\\theta(t)$')
plt.legend();
plt.subplot(2,1,2)
plt.plot(tt,XXC,label="couplée")
plt.plot(tt,sols[:,1],label="simple")
plt.plot(tt,sola[:,1],label="asservie")
plt.plot(tt,XCX,'--',label="optimale")
plt.ylabel('Position $x_c(t)$')
plt.legend();

# tracer du mouvement cde simple
pas = 5
chariots = Chariot(L0,sols[::pas,1],sols[::pas,0],tt[pas])
chariots.trace(titre="commande simple",pas=4)

# animation
anims = chariots.calculanim(-0.5, XF+0.5, -1.2*L0, 0.5)

plt.rc('animation', html='jshtml')
anims
# tracer du mouvement asservisement
chariota = Chariot(L0,sola[::pas,1],sola[::pas,0],tt[pas])
chariota.trace(titre="asservissement simple",pas=4)

# animation
anima = chariota.calculanim(-0.5, XF+0.5, -1.2*L0, 0.5)

plt.rc('animation', html='jshtml')
anima
# tracer du mouvement cde lente couplee
chariotc = Chariot(L0,XXC[::pas],THETA[::pas],tt[pas])
chariotc.trace(titre="cde lente couplée",pas=4)

# animation
animc = chariotc.calculanim(-0.5, XF+0.5, -1.2*L0, 0.5)

plt.rc('animation', html='jshtml')
animc
6.4. Asservissement: regulateur bas niveau#
pour rendre la commande précédente robuste on ajoute un asservissement pour suivre la trajectoire idéale
trajectoire imposée \(x_i(t)\) à partir de \(\xi(t)\) $\(x_i(t) = \xi + \frac{l_0}{g} \ddot{\xi}\)$
commande asservissement (ajout à la commande idéale précédente) $\( F_c = F_i -K_c ( (x_c-x_i) + T_c(\dot{x_c}-\dot{x_i})) \)$
smb = sp.simplify(LM.rhs().subs({g:omega0**2*l0}))
display(smb)
display(xcx)
Kc, Tc = sp.symbols('K_c T_c')
Fca = -Kc*((xc-xcx)+(xcp-xcx.diff(t,1))*Tc)
display(Fca)
smbF = sp.lambdify([theta,xc,thetap,xcp,l0,omega0,alpha,Fc],smb)
FCa = sp.lambdify([t,xc,xcp,tf,xf,omega0,Kc,Tc],Fca)
6.4.1. Simulation#
# parametres
def F(Y,t):
global TF,XF,L0,OMEGA0,ALPHA,KC,TC
FC = 0
# cde idéale
FC += FX(t,XF,TF,OMEGA0,ALPHA)
# asservissement
FC += FCa(t,Y[1],Y[3],TF,XF,OMEGA0,KC,TC)
# 2nd membre
FF = smbF(Y[0],Y[1],Y[2],Y[3],L0,OMEGA0,ALPHA,FC)
return FF[:,0]
# parametres
KC = 200
TC = np.sqrt(4/KC)
print("parametres:",KC,TC,TF)
# simulation
Y0 = [0.0, 0.0, 0., 0. ]
# perturbation
Y0 = [0.01,0.01,0.,0.]
sol = odeint(F,Y0,tt)
THETA = sol[:,0]
XXC = sol[:,1]
print("Erreur sur la position finale")
print("commande couplée erreur xc={:.3f} \t theta={:.3f} deg".format(sol[-1,1]-XF ,np.degrees(sol[-1,0])))
parametres: 200 0.1414213562373095 4.0
Erreur sur la position finale
commande couplée erreur xc=0.002 theta=0.849 deg
plt.subplot(2,1,1)
plt.plot(tt,XXC)
plt.plot(tt,XCX,'--')
plt.ylabel("xc(t)")
plt.subplot(2,1,2)
plt.plot(tt,THETA)
plt.plot(tt,TTHX,'--')
plt.ylabel("$\\theta(t)$")
plt.xlabel('t');

chariot = Chariot(L0,XXC[::pas],THETA[::pas],tt[pas])
chariot.trace(titre="asservisement à commande lente",pas=4)

anim = chariot.calculanim(-0.5, XF+0.5, -1.2*L0, 0.5)

plt.rc('animation', html='jshtml')
anim
6.5. cas avec une longueur différence#
6.5.1. mouvement de la charge#
coordonnees \(\xi(t)\), \(\eta(t)\) dans le repere fixe (voir cours) $\(\xi = l \sin\theta + x-c \)\( \)\(\eta = -l \cos\theta \)$
équations d’équilibre $\( m\ddot\xi = -T \sin\theta\)\( \)\( m\ddot\xi = -T \sin\theta\)$
d’où l’expression de \(\theta\), \(l\) et \(x_c\) en fonction de \(\xi\) et \(\eta\)
trajectoire idéale
on planifie une trajectoire idéale pour \(\xi\) et \(\eta\) , et on en déduit la valeur de \(x_c(t)\) \(\theta(t)\) et \(l(t)\), et le contrôle assoicié
xi, eta = dynamicsymbols('xi eta')
display(sp.Eq(sp.tan(theta), xi.diff(t,t)/(eta.diff(t,t)-g)))
display(sp.Eq(xc, xi - eta*xi.diff(t,t)/(eta.diff(t,t)-g)) )
6.5.1.1. trajectoire planifiée#
trajectoire \(\xi(t)\) tq. $\(\xi(0)=0, \frac{d\xi}{dt}(0)=0, \frac{d^2\xi}{dt^2}(0)=0, \frac{d^3\xi}{dt^3}(0)=0, \frac{d^4\xi}{dt^4}(0)=0\)\( \)\(\xi(t_f)=x_f, \frac{d\xi}{dt}(t_f)=0, \frac{d^2\xi}{dt^2}(t_f)=0, \frac{d^3\xi}{dt^3}(t_f)=0, \frac{d^4\xi}{dt^4}(t_f)=0\)$
\(\xi(t)\) polynome de degré 9 en t
trajectoire \(\eta(t)\) (même condition que \(\xi\)) sauf sur la position
\(\eta=-l_0\) longueur à \(t=0\)
\(\eta=-l_f\) longueur à \(t=t_f\)
xf, tf = sp.symbols('x_f t_f')
x = xf*(t/tf)**5*(126-420*(t/tf)+540*(t/tf)**2-315*(t/tf)**3+70*(t/tf)**4)
display(sp.Eq(xi,x))
# verification
print("dérivées en tf:",x.diff(t,1).subs({t:tf}),x.diff(t,2).subs({t:tf}),x.diff(t,3).subs({t:tf}),x.diff(t,4).subs({t:tf}))
dérivées en tf: 0 0 0 0
lf = sp.symbols('l_f')
y = - l0 - (lf - l0)*(x/xf)
display(sp.Eq(eta,y))
6.5.2. cas linéaire#
on suppose que \(\theta \approx 0\)
trajectoire idéale du chariot et du cable : $\(x_{ref}(t) \approx \xi(t) \)\( \)\(l_{ref}(t) \approx - \eta(t)\)$
d’ou l’asservissement de \(F_c\) et \(C_t\)
xref = x
lref = -y
display(xref, lref)
# asservissement
Kc, Tc = sp.symbols('K_c T_c')
Fca = -Kc*((xc-xref)+(xcp-xref.diff(t,1))*Tc)
Fca = Fca.simplify()
display(Fca)
Kt, Tt = sp.symbols('K_t T_t')
Cta = -g -Kt*((l-lref)+(lp-lref.diff(t,1))*Tt)
Cta = Cta.simplify()
display(Cta)
# fct python
FCa = sp.lambdify([t,xc,xcp,xf,tf,Kc,Tc],Fca,'numpy')
CTa = sp.lambdify([t,l,lp,l0,lf,tf,Kt,Tt,g],Cta,'numpy')
Xref= sp.lambdify([t,xf,tf],xref,'numpy')
Lref= sp.lambdify([t,l0,lf,tf],lref,'numpy')
6.5.3. Equation lagrange avec couple#
alpha = sp.Symbol('alpha')
La = sp.simplify(La.subs({M:alpha*m}))
display(La)
Fc, Ct = sp.symbols('F_c C_t')
LM = LagrangesMethod(La,[theta,xc,l],forcelist=[(C,Fc*alpha*m*RO.x),(Rt,Ct*m*rho*RO.z)], frame = RO)
eq = sp.simplify(LM.form_lagranges_equations())
display(eq)
smb = sp.simplify(sp.simplify(LM.rhs()).doit())
display(smb)
smbF = sp.lambdify([theta,xc,l,thetap,xcp,lp,Fc,Ct,g,alpha],smb,'numpy')
6.5.4. simulation#
# PARAMETRES NUMERIQUES
TF = 10.0
#TF = 8.0
XF = 4.0
ALPHA = 0.5
L0 = 3.0
LF = 2.0
GG = 9.81
KC = 50.0
TC = 20./KC
KT = 50.0
TT = 20./KT
Nt = 400
tt = np.linspace(0,TF,Nt)
# trajectoire idéale
XREF = Xref(tt,XF,TF)
LREF = Lref(tt,L0,LF,TF)
plt.figure(figsize=(12,6))
plt.subplot(1,2,1)
plt.plot(tt,XREF)
plt.title("Position ")
plt.subplot(1,2,2)
plt.plot(tt,LREF)
plt.title("Longueur")
Text(0.5, 1.0, 'Longueur')

# simulation mouvement
def F(Y,t):
global ALPHA,TF,XF
FC = FCa(t,Y[1],Y[4],XF,TF,KC,TC)
CT = CTa(t,Y[2],Y[5],L0,LF,TF,KT,TT,GG)
FF = smbF(Y[0],Y[1],Y[2],Y[3],Y[4],Y[5],FC,CT,GG,ALPHA)
return FF[:,0]
Y0 = [0.0, 0.0, L0, 0., 0., 0. ]
sol = odeint(F,Y0,tt)
THETA = sol[:,0]
XC = sol[:,1]
LL = sol[:,2]
plt.figure(figsize=(12,6))
plt.subplot(1,3,1)
plt.plot(tt,THETA*180/np.pi)
plt.ylim(-10,10)
plt.title("Angle$\\theta$ en degré")
plt.subplot(1,3,2)
plt.plot(tt,XC)
plt.plot(tt,XREF,'--')
plt.title('Position $x_c(t)$')
plt.subplot(1,3,3)
plt.plot(tt,LL)
plt.plot(tt,LREF,'--')
plt.title('longueur $l(t)$')
Text(0.5, 1.0, 'longueur $l(t)$')

# tracer du mouvement
plt.figure(figsize=(12,8))
for i in range(0,len(tt),20):
X,Y = [XC[i],XC[i]+LL[i]*np.sin(THETA[i])],[0.,-LL[i]*np.cos(THETA[i])]
plt.plot(X,Y,'-o',ms=20)
plt.axis('equal')

6.5.5. cas non linéaire#
ddl fonction de la trajectoire (idéale) \(\xi\) \(\eta\) $\( x_{ref} = \xi - \frac{\eta\ddot{\xi}}{\ddot{\eta}-g}\)$
xref = x - y*x.diff(t,2)/(y.diff(t,2)-g)
xref = xref.doit()
lref = sp.sqrt((y*x.diff(t,2)/(y.diff(t,2)-g))**2 + y**2)
lref = lref.doit()
print(xref, lref)
t**5*x_f*(70*t**4/t_f**4 - 315*t**3/t_f**3 + 540*t**2/t_f**2 - 420*t/t_f + 126)/t_f**5 - 10*t**3*x_f*(-l_0 - t**5*(-l_0 + l_f)*(70*t**4/t_f**4 - 315*t**3/t_f**3 + 540*t**2/t_f**2 - 420*t/t_f + 126)/t_f**5)*(140*t**4/t_f**4 - 630*t**3/t_f**3 + 3*t**2*(28*t**2/t_f**2 - 63*t/t_f + 36)/t_f**2 + 1080*t**2/t_f**2 + 5*t*(56*t**3/t_f**3 - 189*t**2/t_f**2 + 216*t/t_f - 84)/t_f - 840*t/t_f + 252)/(t_f**5*(-g + 10*t**3*(l_0 - l_f)*(140*t**4/t_f**4 - 630*t**3/t_f**3 + 3*t**2*(28*t**2/t_f**2 - 63*t/t_f + 36)/t_f**2 + 1080*t**2/t_f**2 + 5*t*(56*t**3/t_f**3 - 189*t**2/t_f**2 + 216*t/t_f - 84)/t_f - 840*t/t_f + 252)/t_f**5)) sqrt(100*t**6*x_f**2*(-l_0 - t**5*(-l_0 + l_f)*(70*t**4/t_f**4 - 315*t**3/t_f**3 + 540*t**2/t_f**2 - 420*t/t_f + 126)/t_f**5)**2*(140*t**4/t_f**4 - 630*t**3/t_f**3 + 3*t**2*(28*t**2/t_f**2 - 63*t/t_f + 36)/t_f**2 + 1080*t**2/t_f**2 + 5*t*(56*t**3/t_f**3 - 189*t**2/t_f**2 + 216*t/t_f - 84)/t_f - 840*t/t_f + 252)**2/(t_f**10*(-g + 10*t**3*(l_0 - l_f)*(140*t**4/t_f**4 - 630*t**3/t_f**3 + 3*t**2*(28*t**2/t_f**2 - 63*t/t_f + 36)/t_f**2 + 1080*t**2/t_f**2 + 5*t*(56*t**3/t_f**3 - 189*t**2/t_f**2 + 216*t/t_f - 84)/t_f - 840*t/t_f + 252)/t_f**5)**2) + (-l_0 - t**5*(-l_0 + l_f)*(70*t**4/t_f**4 - 315*t**3/t_f**3 + 540*t**2/t_f**2 - 420*t/t_f + 126)/t_f**5)**2)
# asservissement pour obtenir la trajectoire
Kc, Tc = sp.symbols('K_c T_c')
Fca = -Kc*((xc-xref)+(xcp-xref.diff(t,1))*Tc)
#Fca = Fca.simplify()
#display(Fca)
print(Fca)
Kt, Tt = sp.symbols('K_t T_t')
Cta = -g - Kt*((l-lref)+(lp-lref.diff(t,1))*Tt)
#Cta = Cta.simplify()
#display(Cta)
print(Cta)
-K_c*(T_c*(-t**5*x_f*(280*t**3/t_f**4 - 945*t**2/t_f**3 + 1080*t/t_f**2 - 420/t_f)/t_f**5 - 5*t**4*x_f*(70*t**4/t_f**4 - 315*t**3/t_f**3 + 540*t**2/t_f**2 - 420*t/t_f + 126)/t_f**5 + 10*t**3*x_f*(-l_0 - t**5*(-l_0 + l_f)*(70*t**4/t_f**4 - 315*t**3/t_f**3 + 540*t**2/t_f**2 - 420*t/t_f + 126)/t_f**5)*(560*t**3/t_f**4 + 3*t**2*(56*t/t_f**2 - 63/t_f)/t_f**2 - 1890*t**2/t_f**3 + 5*t*(168*t**2/t_f**3 - 378*t/t_f**2 + 216/t_f)/t_f + 6*t*(28*t**2/t_f**2 - 63*t/t_f + 36)/t_f**2 + 2160*t/t_f**2 + 5*(56*t**3/t_f**3 - 189*t**2/t_f**2 + 216*t/t_f - 84)/t_f - 840/t_f)/(t_f**5*(-g + 10*t**3*(l_0 - l_f)*(140*t**4/t_f**4 - 630*t**3/t_f**3 + 3*t**2*(28*t**2/t_f**2 - 63*t/t_f + 36)/t_f**2 + 1080*t**2/t_f**2 + 5*t*(56*t**3/t_f**3 - 189*t**2/t_f**2 + 216*t/t_f - 84)/t_f - 840*t/t_f + 252)/t_f**5)) + 10*t**3*x_f*(-t**5*(-l_0 + l_f)*(280*t**3/t_f**4 - 945*t**2/t_f**3 + 1080*t/t_f**2 - 420/t_f)/t_f**5 - 5*t**4*(-l_0 + l_f)*(70*t**4/t_f**4 - 315*t**3/t_f**3 + 540*t**2/t_f**2 - 420*t/t_f + 126)/t_f**5)*(140*t**4/t_f**4 - 630*t**3/t_f**3 + 3*t**2*(28*t**2/t_f**2 - 63*t/t_f + 36)/t_f**2 + 1080*t**2/t_f**2 + 5*t*(56*t**3/t_f**3 - 189*t**2/t_f**2 + 216*t/t_f - 84)/t_f - 840*t/t_f + 252)/(t_f**5*(-g + 10*t**3*(l_0 - l_f)*(140*t**4/t_f**4 - 630*t**3/t_f**3 + 3*t**2*(28*t**2/t_f**2 - 63*t/t_f + 36)/t_f**2 + 1080*t**2/t_f**2 + 5*t*(56*t**3/t_f**3 - 189*t**2/t_f**2 + 216*t/t_f - 84)/t_f - 840*t/t_f + 252)/t_f**5)) + 10*t**3*x_f*(-l_0 - t**5*(-l_0 + l_f)*(70*t**4/t_f**4 - 315*t**3/t_f**3 + 540*t**2/t_f**2 - 420*t/t_f + 126)/t_f**5)*(-10*t**3*(l_0 - l_f)*(560*t**3/t_f**4 + 3*t**2*(56*t/t_f**2 - 63/t_f)/t_f**2 - 1890*t**2/t_f**3 + 5*t*(168*t**2/t_f**3 - 378*t/t_f**2 + 216/t_f)/t_f + 6*t*(28*t**2/t_f**2 - 63*t/t_f + 36)/t_f**2 + 2160*t/t_f**2 + 5*(56*t**3/t_f**3 - 189*t**2/t_f**2 + 216*t/t_f - 84)/t_f - 840/t_f)/t_f**5 - 30*t**2*(l_0 - l_f)*(140*t**4/t_f**4 - 630*t**3/t_f**3 + 3*t**2*(28*t**2/t_f**2 - 63*t/t_f + 36)/t_f**2 + 1080*t**2/t_f**2 + 5*t*(56*t**3/t_f**3 - 189*t**2/t_f**2 + 216*t/t_f - 84)/t_f - 840*t/t_f + 252)/t_f**5)*(140*t**4/t_f**4 - 630*t**3/t_f**3 + 3*t**2*(28*t**2/t_f**2 - 63*t/t_f + 36)/t_f**2 + 1080*t**2/t_f**2 + 5*t*(56*t**3/t_f**3 - 189*t**2/t_f**2 + 216*t/t_f - 84)/t_f - 840*t/t_f + 252)/(t_f**5*(-g + 10*t**3*(l_0 - l_f)*(140*t**4/t_f**4 - 630*t**3/t_f**3 + 3*t**2*(28*t**2/t_f**2 - 63*t/t_f + 36)/t_f**2 + 1080*t**2/t_f**2 + 5*t*(56*t**3/t_f**3 - 189*t**2/t_f**2 + 216*t/t_f - 84)/t_f - 840*t/t_f + 252)/t_f**5)**2) + 30*t**2*x_f*(-l_0 - t**5*(-l_0 + l_f)*(70*t**4/t_f**4 - 315*t**3/t_f**3 + 540*t**2/t_f**2 - 420*t/t_f + 126)/t_f**5)*(140*t**4/t_f**4 - 630*t**3/t_f**3 + 3*t**2*(28*t**2/t_f**2 - 63*t/t_f + 36)/t_f**2 + 1080*t**2/t_f**2 + 5*t*(56*t**3/t_f**3 - 189*t**2/t_f**2 + 216*t/t_f - 84)/t_f - 840*t/t_f + 252)/(t_f**5*(-g + 10*t**3*(l_0 - l_f)*(140*t**4/t_f**4 - 630*t**3/t_f**3 + 3*t**2*(28*t**2/t_f**2 - 63*t/t_f + 36)/t_f**2 + 1080*t**2/t_f**2 + 5*t*(56*t**3/t_f**3 - 189*t**2/t_f**2 + 216*t/t_f - 84)/t_f - 840*t/t_f + 252)/t_f**5)) + Derivative(x_c(t), t)) - t**5*x_f*(70*t**4/t_f**4 - 315*t**3/t_f**3 + 540*t**2/t_f**2 - 420*t/t_f + 126)/t_f**5 + 10*t**3*x_f*(-l_0 - t**5*(-l_0 + l_f)*(70*t**4/t_f**4 - 315*t**3/t_f**3 + 540*t**2/t_f**2 - 420*t/t_f + 126)/t_f**5)*(140*t**4/t_f**4 - 630*t**3/t_f**3 + 3*t**2*(28*t**2/t_f**2 - 63*t/t_f + 36)/t_f**2 + 1080*t**2/t_f**2 + 5*t*(56*t**3/t_f**3 - 189*t**2/t_f**2 + 216*t/t_f - 84)/t_f - 840*t/t_f + 252)/(t_f**5*(-g + 10*t**3*(l_0 - l_f)*(140*t**4/t_f**4 - 630*t**3/t_f**3 + 3*t**2*(28*t**2/t_f**2 - 63*t/t_f + 36)/t_f**2 + 1080*t**2/t_f**2 + 5*t*(56*t**3/t_f**3 - 189*t**2/t_f**2 + 216*t/t_f - 84)/t_f - 840*t/t_f + 252)/t_f**5)) + x_c(t))
-K_t*(T_t*(Derivative(l(t), t) - (50*t**6*x_f**2*(-l_0 - t**5*(-l_0 + l_f)*(70*t**4/t_f**4 - 315*t**3/t_f**3 + 540*t**2/t_f**2 - 420*t/t_f + 126)/t_f**5)**2*(140*t**4/t_f**4 - 630*t**3/t_f**3 + 3*t**2*(28*t**2/t_f**2 - 63*t/t_f + 36)/t_f**2 + 1080*t**2/t_f**2 + 5*t*(56*t**3/t_f**3 - 189*t**2/t_f**2 + 216*t/t_f - 84)/t_f - 840*t/t_f + 252)*(1120*t**3/t_f**4 + 6*t**2*(56*t/t_f**2 - 63/t_f)/t_f**2 - 3780*t**2/t_f**3 + 10*t*(168*t**2/t_f**3 - 378*t/t_f**2 + 216/t_f)/t_f + 12*t*(28*t**2/t_f**2 - 63*t/t_f + 36)/t_f**2 + 4320*t/t_f**2 + 10*(56*t**3/t_f**3 - 189*t**2/t_f**2 + 216*t/t_f - 84)/t_f - 1680/t_f)/(t_f**10*(-g + 10*t**3*(l_0 - l_f)*(140*t**4/t_f**4 - 630*t**3/t_f**3 + 3*t**2*(28*t**2/t_f**2 - 63*t/t_f + 36)/t_f**2 + 1080*t**2/t_f**2 + 5*t*(56*t**3/t_f**3 - 189*t**2/t_f**2 + 216*t/t_f - 84)/t_f - 840*t/t_f + 252)/t_f**5)**2) + 50*t**6*x_f**2*(-l_0 - t**5*(-l_0 + l_f)*(70*t**4/t_f**4 - 315*t**3/t_f**3 + 540*t**2/t_f**2 - 420*t/t_f + 126)/t_f**5)*(-2*t**5*(-l_0 + l_f)*(280*t**3/t_f**4 - 945*t**2/t_f**3 + 1080*t/t_f**2 - 420/t_f)/t_f**5 - 10*t**4*(-l_0 + l_f)*(70*t**4/t_f**4 - 315*t**3/t_f**3 + 540*t**2/t_f**2 - 420*t/t_f + 126)/t_f**5)*(140*t**4/t_f**4 - 630*t**3/t_f**3 + 3*t**2*(28*t**2/t_f**2 - 63*t/t_f + 36)/t_f**2 + 1080*t**2/t_f**2 + 5*t*(56*t**3/t_f**3 - 189*t**2/t_f**2 + 216*t/t_f - 84)/t_f - 840*t/t_f + 252)**2/(t_f**10*(-g + 10*t**3*(l_0 - l_f)*(140*t**4/t_f**4 - 630*t**3/t_f**3 + 3*t**2*(28*t**2/t_f**2 - 63*t/t_f + 36)/t_f**2 + 1080*t**2/t_f**2 + 5*t*(56*t**3/t_f**3 - 189*t**2/t_f**2 + 216*t/t_f - 84)/t_f - 840*t/t_f + 252)/t_f**5)**2) + 50*t**6*x_f**2*(-l_0 - t**5*(-l_0 + l_f)*(70*t**4/t_f**4 - 315*t**3/t_f**3 + 540*t**2/t_f**2 - 420*t/t_f + 126)/t_f**5)**2*(-20*t**3*(l_0 - l_f)*(560*t**3/t_f**4 + 3*t**2*(56*t/t_f**2 - 63/t_f)/t_f**2 - 1890*t**2/t_f**3 + 5*t*(168*t**2/t_f**3 - 378*t/t_f**2 + 216/t_f)/t_f + 6*t*(28*t**2/t_f**2 - 63*t/t_f + 36)/t_f**2 + 2160*t/t_f**2 + 5*(56*t**3/t_f**3 - 189*t**2/t_f**2 + 216*t/t_f - 84)/t_f - 840/t_f)/t_f**5 - 60*t**2*(l_0 - l_f)*(140*t**4/t_f**4 - 630*t**3/t_f**3 + 3*t**2*(28*t**2/t_f**2 - 63*t/t_f + 36)/t_f**2 + 1080*t**2/t_f**2 + 5*t*(56*t**3/t_f**3 - 189*t**2/t_f**2 + 216*t/t_f - 84)/t_f - 840*t/t_f + 252)/t_f**5)*(140*t**4/t_f**4 - 630*t**3/t_f**3 + 3*t**2*(28*t**2/t_f**2 - 63*t/t_f + 36)/t_f**2 + 1080*t**2/t_f**2 + 5*t*(56*t**3/t_f**3 - 189*t**2/t_f**2 + 216*t/t_f - 84)/t_f - 840*t/t_f + 252)**2/(t_f**10*(-g + 10*t**3*(l_0 - l_f)*(140*t**4/t_f**4 - 630*t**3/t_f**3 + 3*t**2*(28*t**2/t_f**2 - 63*t/t_f + 36)/t_f**2 + 1080*t**2/t_f**2 + 5*t*(56*t**3/t_f**3 - 189*t**2/t_f**2 + 216*t/t_f - 84)/t_f - 840*t/t_f + 252)/t_f**5)**3) + 300*t**5*x_f**2*(-l_0 - t**5*(-l_0 + l_f)*(70*t**4/t_f**4 - 315*t**3/t_f**3 + 540*t**2/t_f**2 - 420*t/t_f + 126)/t_f**5)**2*(140*t**4/t_f**4 - 630*t**3/t_f**3 + 3*t**2*(28*t**2/t_f**2 - 63*t/t_f + 36)/t_f**2 + 1080*t**2/t_f**2 + 5*t*(56*t**3/t_f**3 - 189*t**2/t_f**2 + 216*t/t_f - 84)/t_f - 840*t/t_f + 252)**2/(t_f**10*(-g + 10*t**3*(l_0 - l_f)*(140*t**4/t_f**4 - 630*t**3/t_f**3 + 3*t**2*(28*t**2/t_f**2 - 63*t/t_f + 36)/t_f**2 + 1080*t**2/t_f**2 + 5*t*(56*t**3/t_f**3 - 189*t**2/t_f**2 + 216*t/t_f - 84)/t_f - 840*t/t_f + 252)/t_f**5)**2) + (-l_0 - t**5*(-l_0 + l_f)*(70*t**4/t_f**4 - 315*t**3/t_f**3 + 540*t**2/t_f**2 - 420*t/t_f + 126)/t_f**5)*(-2*t**5*(-l_0 + l_f)*(280*t**3/t_f**4 - 945*t**2/t_f**3 + 1080*t/t_f**2 - 420/t_f)/t_f**5 - 10*t**4*(-l_0 + l_f)*(70*t**4/t_f**4 - 315*t**3/t_f**3 + 540*t**2/t_f**2 - 420*t/t_f + 126)/t_f**5)/2)/sqrt(100*t**6*x_f**2*(-l_0 - t**5*(-l_0 + l_f)*(70*t**4/t_f**4 - 315*t**3/t_f**3 + 540*t**2/t_f**2 - 420*t/t_f + 126)/t_f**5)**2*(140*t**4/t_f**4 - 630*t**3/t_f**3 + 3*t**2*(28*t**2/t_f**2 - 63*t/t_f + 36)/t_f**2 + 1080*t**2/t_f**2 + 5*t*(56*t**3/t_f**3 - 189*t**2/t_f**2 + 216*t/t_f - 84)/t_f - 840*t/t_f + 252)**2/(t_f**10*(-g + 10*t**3*(l_0 - l_f)*(140*t**4/t_f**4 - 630*t**3/t_f**3 + 3*t**2*(28*t**2/t_f**2 - 63*t/t_f + 36)/t_f**2 + 1080*t**2/t_f**2 + 5*t*(56*t**3/t_f**3 - 189*t**2/t_f**2 + 216*t/t_f - 84)/t_f - 840*t/t_f + 252)/t_f**5)**2) + (-l_0 - t**5*(-l_0 + l_f)*(70*t**4/t_f**4 - 315*t**3/t_f**3 + 540*t**2/t_f**2 - 420*t/t_f + 126)/t_f**5)**2)) - sqrt(100*t**6*x_f**2*(-l_0 - t**5*(-l_0 + l_f)*(70*t**4/t_f**4 - 315*t**3/t_f**3 + 540*t**2/t_f**2 - 420*t/t_f + 126)/t_f**5)**2*(140*t**4/t_f**4 - 630*t**3/t_f**3 + 3*t**2*(28*t**2/t_f**2 - 63*t/t_f + 36)/t_f**2 + 1080*t**2/t_f**2 + 5*t*(56*t**3/t_f**3 - 189*t**2/t_f**2 + 216*t/t_f - 84)/t_f - 840*t/t_f + 252)**2/(t_f**10*(-g + 10*t**3*(l_0 - l_f)*(140*t**4/t_f**4 - 630*t**3/t_f**3 + 3*t**2*(28*t**2/t_f**2 - 63*t/t_f + 36)/t_f**2 + 1080*t**2/t_f**2 + 5*t*(56*t**3/t_f**3 - 189*t**2/t_f**2 + 216*t/t_f - 84)/t_f - 840*t/t_f + 252)/t_f**5)**2) + (-l_0 - t**5*(-l_0 + l_f)*(70*t**4/t_f**4 - 315*t**3/t_f**3 + 540*t**2/t_f**2 - 420*t/t_f + 126)/t_f**5)**2) + l(t)) - g
# fct numpy
FCa = sp.lambdify([t,xc,xcp,xf,l0,lf,tf,Kc,Tc,g],Fca,'numpy')
CTa = sp.lambdify([t,l,lp,xf,l0,lf,tf,Kt,Tt,g],Cta,'numpy')
Xref= sp.lambdify([t,xf,l0,lf,tf,g],xref,'numpy')
Lref= sp.lambdify([t,xf,l0,lf,tf,g],lref,'numpy')
# trajectoire ideale
XREF = Xref(tt,XF,L0,LF,TF,GG)
LREF = Lref(tt,XF,L0,LF,TF,GG)
plt.figure(figsize=(12,6))
plt.subplot(1,2,1)
plt.plot(tt,XREF)
plt.title("Position ")
plt.subplot(1,2,2)
plt.plot(tt,LREF)
plt.title("Longueur")
Text(0.5, 1.0, 'Longueur')

# simulation mouvement
def F(Y,t):
global ALPHA,TF,XF
FC = FCa(t,Y[1],Y[4],XF,L0,LF,TF,KC,TC,GG)
CT = CTa(t,Y[2],Y[5],XF,L0,LF,TF,KT,TT,GG)
FF = smbF(Y[0],Y[1],Y[2],Y[3],Y[4],Y[5],FC,CT,GG,ALPHA)
return FF[:,0]
Y0 = [0.0, 0.0, L0, 0., 0., 0. ]
print(F(Y0,0.))
sol = odeint(F,Y0,tt)
THETA = sol[:,0]
XC = sol[:,1]
LL = sol[:,2]
[0. 0. 0. 0. 0. 0.]
plt.figure(figsize=(12,6))
plt.subplot(1,3,1)
plt.plot(tt,THETA*180/np.pi)
plt.ylim(-10,10)
plt.title("Angle$\\theta$ en degré")
plt.subplot(1,3,2)
plt.plot(tt,XC)
plt.plot(tt,XREF,'--')
plt.title('Position $x_c(t)$')
plt.subplot(1,3,3)
plt.plot(tt,LL)
plt.plot(tt,LREF,'--')
plt.title('longueur $l(t)$')
Text(0.5, 1.0, 'longueur $l(t)$')

# tracer du mouvement
plt.figure(figsize=(12,8))
for i in range(0,len(tt),20):
X,Y = [XC[i],XC[i]+LL[i]*np.sin(THETA[i])],[0.,-LL[i]*np.cos(THETA[i])]
plt.plot(X,Y,'-o',ms=20)
plt.axis('equal')

6.6. Evitement d’un obstacle#
au point x=xref/2 on passe au dessus d’un obstacle pour avoir une longueur la < min(l0,lf)
sans obstacle la trajectoire droite: $\( y = l_0 + (l_f-l_0) x \)$
avec obstacle tq \(y=y_a\): on xie par une parabole $\( y = l_0 + (l_f-l_0) x (a+bx+cx^2)\)$
avec 3 conditions:
\(y = l_f\) en x=1
\(y = l_a\) en x=1/2
\(\dot{y} = 0\) en x=1/2
pour \(l_a = 2l_f - l_0\) on trouve: