1515
1616SHOW_ANIMATION = True
1717
18- MAX_TIME = 100.0 # Maximum simulation time
19- DT = 0.1 # Time tick
2018
19+ class LQRPlanner :
2120
22- def LQRplanning (sx , sy , gx , gy , show_animation = SHOW_ANIMATION ):
21+ def __init__ (self ):
22+ self .MAX_TIME = 100.0 # Maximum simulation time
23+ self .DT = 0.1 # Time tick
24+ self .GOAL_DIST = 0.1
25+ self .MAX_ITER = 150
26+ self .EPS = 0.01
2327
24- rx , ry = [ sx ], [ sy ]
28+ def lqr_planning ( self , sx , sy , gx , gy , show_animation = SHOW_ANIMATION ):
2529
26- x = np . array ([ sx - gx , sy - gy ]). reshape ( 2 , 1 ) # State vector
30+ rx , ry = [ sx ], [ sy ]
2731
28- # Linear system model
29- A , B = get_system_model ()
32+ x = np .array ([sx - gx , sy - gy ]).reshape (2 , 1 ) # State vector
3033
31- found_path = False
34+ # Linear system model
35+ A , B = self .get_system_model ()
3236
33- time = 0.0
34- while time <= MAX_TIME :
35- time += DT
37+ found_path = False
3638
37- u = LQR_control (A , B , x )
39+ time = 0.0
40+ while time <= self .MAX_TIME :
41+ time += self .DT
3842
39- x = A @ x + B @ u
43+ u = self . lqr_control ( A , B , x )
4044
41- rx .append (x [0 , 0 ] + gx )
42- ry .append (x [1 , 0 ] + gy )
45+ x = A @ x + B @ u
4346
44- d = math .sqrt ((gx - rx [- 1 ])** 2 + (gy - ry [- 1 ])** 2 )
45- if d <= 0.1 :
46- # print("Goal!!")
47- found_path = True
48- break
47+ rx .append (x [0 , 0 ] + gx )
48+ ry .append (x [1 , 0 ] + gy )
4949
50- # animation
51- if show_animation : # pragma: no cover
52- plt .plot (sx , sy , "or" )
53- plt .plot (gx , gy , "ob" )
54- plt .plot (rx , ry , "-r" )
55- plt .axis ("equal" )
56- plt .pause (1.0 )
50+ d = math .sqrt ((gx - rx [- 1 ]) ** 2 + (gy - ry [- 1 ]) ** 2 )
51+ if d <= self .GOAL_DIST :
52+ found_path = True
53+ break
5754
58- if not found_path :
59- print ("Cannot found path" )
60- return [], []
55+ # animation
56+ if show_animation : # pragma: no cover
57+ plt .plot (sx , sy , "or" )
58+ plt .plot (gx , gy , "ob" )
59+ plt .plot (rx , ry , "-r" )
60+ plt .axis ("equal" )
61+ plt .pause (1.0 )
6162
62- return rx , ry
63+ if not found_path :
64+ print ("Cannot found path" )
65+ return [], []
6366
67+ return rx , ry
6468
65- def solve_DARE (A , B , Q , R ):
66- """
67- solve a discrete time_Algebraic Riccati equation (DARE)
68- """
69- X = Q
70- maxiter = 150
71- eps = 0.01
69+ def solve_dare (self , A , B , Q , R ):
70+ """
71+ solve a discrete time_Algebraic Riccati equation (DARE)
72+ """
73+ X , Xn = Q , Q
7274
73- for i in range (maxiter ):
74- Xn = A .T * X * A - A .T * X * B * \
75- la .inv (R + B .T * X * B ) * B .T * X * A + Q
76- if (abs (Xn - X )).max () < eps :
77- break
78- X = Xn
75+ for i in range (self . MAX_ITER ):
76+ Xn = A .T * X * A - A .T * X * B * \
77+ la .inv (R + B .T * X * B ) * B .T * X * A + Q
78+ if (abs (Xn - X )).max () < self . EPS :
79+ break
80+ X = Xn
7981
80- return Xn
82+ return Xn
8183
84+ def dlqr (self , A , B , Q , R ):
85+ """Solve the discrete time lqr controller.
86+ x[k+1] = A x[k] + B u[k]
87+ cost = sum x[k].T*Q*x[k] + u[k].T*R*u[k]
88+ # ref Bertsekas, p.151
89+ """
8290
83- def dlqr (A , B , Q , R ):
84- """Solve the discrete time lqr controller.
85- x[k+1] = A x[k] + B u[k]
86- cost = sum x[k].T*Q*x[k] + u[k].T*R*u[k]
87- # ref Bertsekas, p.151
88- """
91+ # first, try to solve the ricatti equation
92+ X = self .solve_dare (A , B , Q , R )
8993
90- # first, try to solve the ricatti equation
91- X = solve_DARE ( A , B , Q , R )
94+ # compute the LQR gain
95+ K = la . inv ( B . T @ X @ B + R ) @ ( B . T @ X @ A )
9296
93- # compute the LQR gain
94- K = la .inv (B .T @ X @ B + R ) @ (B .T @ X @ A )
97+ eigValues = la .eigvals (A - B @ K )
9598
96- eigVals , eigVecs = la . eig ( A - B @ K )
99+ return K , X , eigValues
97100
98- return K , X , eigVals
101+ def get_system_model ( self ):
99102
103+ A = np .array ([[self .DT , 1.0 ],
104+ [0.0 , self .DT ]])
105+ B = np .array ([0.0 , 1.0 ]).reshape (2 , 1 )
100106
101- def get_system_model ():
107+ return A , B
102108
103- A = np .array ([[DT , 1.0 ],
104- [0.0 , DT ]])
105- B = np .array ([0.0 , 1.0 ]).reshape (2 , 1 )
109+ def lqr_control (self , A , B , x ):
106110
107- return A , B
111+ Kopt , X , ev = self . dlqr ( A , B , np . eye ( 2 ), np . eye ( 1 ))
108112
113+ u = - Kopt @ x
109114
110- def LQR_control (A , B , x ):
111-
112- Kopt , X , ev = dlqr (A , B , np .eye (2 ), np .eye (1 ))
113-
114- u = - Kopt @ x
115-
116- return u
115+ return u
117116
118117
119118def main ():
@@ -122,13 +121,15 @@ def main():
122121 ntest = 10 # number of goal
123122 area = 100.0 # sampling area
124123
124+ lqr_planner = LQRPlanner ()
125+
125126 for i in range (ntest ):
126127 sx = 6.0
127128 sy = 6.0
128129 gx = random .uniform (- area , area )
129130 gy = random .uniform (- area , area )
130131
131- rx , ry = LQRplanning (sx , sy , gx , gy )
132+ rx , ry = lqr_planner . lqr_planning (sx , sy , gx , gy )
132133
133134 if SHOW_ANIMATION : # pragma: no cover
134135 plt .plot (sx , sy , "or" )
0 commit comments