@@ -43,7 +43,7 @@ class Rocket_Model_6DoF:
4343    A 6 degree of freedom rocket landing problem. 
4444    """ 
4545
46-     def  __init__ (self ):
46+     def  __init__ (self ,  rng ):
4747        """ 
4848        A large r_scale for a small scale problem will 
4949        ead to numerical problems as parameters become excessively small 
@@ -92,7 +92,7 @@ def __init__(self):
9292        # Vector from thrust point to CoM 
9393        self .r_T_B  =  np .array ([- 1e-2 , 0. , 0. ])
9494
95-         self .set_random_initial_state ()
95+         self .set_random_initial_state (rng )
9696
9797        self .x_init  =  np .concatenate (
9898            ((self .m_wet ,), self .r_I_init , self .v_I_init , self .q_B_I_init , self .w_B_init ))
@@ -102,29 +102,32 @@ def __init__(self):
102102        self .r_scale  =  np .linalg .norm (self .r_I_init )
103103        self .m_scale  =  self .m_wet 
104104
105-     def  set_random_initial_state (self ):
105+     def  set_random_initial_state (self , rng ):
106+         if  rng  is  None :
107+             rng  =  np .random .default_rng ()
108+ 
106109        self .r_I_init  =  np .array ((0. , 0. , 0. ))
107-         self .r_I_init [0 ] =  np . random .uniform (3 , 4 )
108-         self .r_I_init [1 :3 ] =  np . random .uniform (- 2 , 2 , size = 2 )
110+         self .r_I_init [0 ] =  rng .uniform (3 , 4 )
111+         self .r_I_init [1 :3 ] =  rng .uniform (- 2 , 2 , size = 2 )
109112
110113        self .v_I_init  =  np .array ((0. , 0. , 0. ))
111-         self .v_I_init [0 ] =  np . random .uniform (- 1 , - 0.5 )
112-         self .v_I_init [1 :3 ] =  np . random . uniform (
113-             - 0.5 ,  - 0.2 ,  size = 2 ) *  self .r_I_init [1 :3 ]
114+         self .v_I_init [0 ] =  rng .uniform (- 1 , - 0.5 )
115+         self .v_I_init [1 :3 ] =  rng . uniform (- 0.5 ,  - 0.2 , 
116+                                           size = 2 ) *  self .r_I_init [1 :3 ]
114117
115118        self .q_B_I_init  =  self .euler_to_quat ((0 ,
116-                                               np . random .uniform (- 30 , 30 ),
117-                                               np . random .uniform (- 30 , 30 )))
119+                                               rng .uniform (- 30 , 30 ),
120+                                               rng .uniform (- 30 , 30 )))
118121        self .w_B_init  =  np .deg2rad ((0 ,
119-                                     np . random .uniform (- 20 , 20 ),
120-                                     np . random .uniform (- 20 , 20 )))
122+                                     rng .uniform (- 20 , 20 ),
123+                                     rng .uniform (- 20 , 20 )))
121124
122125    def  f_func (self , x , u ):
123126        m , rx , ry , rz , vx , vy , vz , q0 , q1 , q2 , q3 , wx , wy , wz  =  x [0 ], x [1 ], x [
124127            2 ], x [3 ], x [4 ], x [5 ], x [6 ], x [7 ], x [8 ], x [9 ], x [10 ], x [11 ], x [12 ], x [13 ]
125128        ux , uy , uz  =  u [0 ], u [1 ], u [2 ]
126129
127-         return  np .matrix ([
130+         return  np .array ([
128131            [- 0.01  *  np .sqrt (ux ** 2  +  uy ** 2  +  uz ** 2 )],
129132            [vx ],
130133            [vy ],
@@ -149,7 +152,7 @@ def A_func(self, x, u):
149152            2 ], x [3 ], x [4 ], x [5 ], x [6 ], x [7 ], x [8 ], x [9 ], x [10 ], x [11 ], x [12 ], x [13 ]
150153        ux , uy , uz  =  u [0 ], u [1 ], u [2 ]
151154
152-         return  np .matrix ([
155+         return  np .array ([
153156            [0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ],
154157            [0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ],
155158            [0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ],
@@ -177,7 +180,7 @@ def B_func(self, x, u):
177180            2 ], x [3 ], x [4 ], x [5 ], x [6 ], x [7 ], x [8 ], x [9 ], x [10 ], x [11 ], x [12 ], x [13 ]
178181        ux , uy , uz  =  u [0 ], u [1 ], u [2 ]
179182
180-         return  np .matrix ([
183+         return  np .array ([
181184            [- 0.01  *  ux  /  np .sqrt (ux ** 2  +  uy ** 2  +  uz ** 2 ),
182185             - 0.01  *  uy  /  np .sqrt (ux  **  2  +  uy ** 2  +  uz ** 2 ),
183186             - 0.01  *  uz  /  np .sqrt (ux ** 2  +  uy ** 2  +  uz ** 2 )],
@@ -219,14 +222,14 @@ def euler_to_quat(self, a):
219222        return  q 
220223
221224    def  skew (self , v ):
222-         return  np .matrix ([
225+         return  np .array ([
223226            [0 , - v [2 ], v [1 ]],
224227            [v [2 ], 0 , - v [0 ]],
225228            [- v [1 ], v [0 ], 0 ]
226229        ])
227230
228231    def  dir_cosine (self , q ):
229-         return  np .matrix ([
232+         return  np .array ([
230233            [1  -  2  *  (q [2 ] **  2  +  q [3 ] **  2 ), 2  *  (q [1 ] *  q [2 ]
231234                                                   +  q [0 ] *  q [3 ]), 2  *  (q [1 ] *  q [3 ] -  q [0 ] *  q [2 ])],
232235            [2  *  (q [1 ] *  q [2 ] -  q [0 ] *  q [3 ]), 1  -  2 
@@ -236,7 +239,7 @@ def dir_cosine(self, q):
236239        ])
237240
238241    def  omega (self , w ):
239-         return  np .matrix ([
242+         return  np .array ([
240243            [0 , - w [0 ], - w [1 ], - w [2 ]],
241244            [w [0 ], 0 , w [2 ], - w [1 ]],
242245            [w [1 ], - w [2 ], 0 , w [0 ]],
@@ -304,7 +307,7 @@ def get_constraints(self, X_v, U_v, X_last_p, U_last_p):
304307        ]
305308
306309        # linearized lower thrust constraint 
307-         rhs  =  [U_last_p [:, k ] /  cvxpy .norm (U_last_p [:, k ]) *  U_v [:, k ]
310+         rhs  =  [U_last_p [:, k ] /  cvxpy .norm (U_last_p [:, k ]) @  U_v [:, k ]
308311               for  k  in  range (X_v .shape [1 ])]
309312        constraints  +=  [
310313            self .T_min  <=  cvxpy .vstack (rhs )
@@ -460,11 +463,11 @@ def __init__(self, m, K):
460463        # x_t+1 = A_*x_t+B_*U_t+C_*U_T+1*S_*sigma+zbar+nu 
461464        constraints  +=  [
462465            self .var ['X' ][:, k  +  1 ] == 
463-             cvxpy .reshape (self .par ['A_bar' ][:, k ], (m .n_x , m .n_x )) * 
466+             cvxpy .reshape (self .par ['A_bar' ][:, k ], (m .n_x , m .n_x )) @ 
464467            self .var ['X' ][:, k ] + 
465-             cvxpy .reshape (self .par ['B_bar' ][:, k ], (m .n_x , m .n_u )) * 
468+             cvxpy .reshape (self .par ['B_bar' ][:, k ], (m .n_x , m .n_u )) @ 
466469            self .var ['U' ][:, k ] + 
467-             cvxpy .reshape (self .par ['C_bar' ][:, k ], (m .n_x , m .n_u )) * 
470+             cvxpy .reshape (self .par ['C_bar' ][:, k ], (m .n_x , m .n_u )) @ 
468471            self .var ['U' ][:, k  +  1 ] + 
469472            self .par ['S_bar' ][:, k ] *  self .var ['sigma' ] + 
470473            self .par ['z_bar' ][:, k ] + 
@@ -606,9 +609,9 @@ def plot_animation(X, U):  # pragma: no cover
606609        plt .pause (0.5 )
607610
608611
609- def  main ():
612+ def  main (rng = None ):
610613    print ("start!!" )
611-     m  =  Rocket_Model_6DoF ()
614+     m  =  Rocket_Model_6DoF (rng )
612615
613616    # state and input list 
614617    X  =  np .empty (shape = [m .n_x , K ])
0 commit comments