|
| 1 | +# # Monte Carlo class usage |
| 2 | + |
| 3 | +# Finally the Monte Carlo simulations can be performed using a dedicated class called `MonteCarlo`. |
| 4 | +# Say goodbye to the long and tedious process of creating the Monte Carlo Simulations throughout jupyter notebooks. |
| 5 | +# |
| 6 | + |
| 7 | +# First, let's import the necessary libraries, including the newest `MonteCarlo` class! |
| 8 | +# |
| 9 | + |
| 10 | +from rocketpy import Environment, SolidMotor, Rocket, Flight, MonteCarlo, Function |
| 11 | +from rocketpy.stochastic import ( |
| 12 | + StochasticEnvironment, |
| 13 | + StochasticSolidMotor, |
| 14 | + StochasticRocket, |
| 15 | + StochasticFlight, |
| 16 | + StochasticNoseCone, |
| 17 | + StochasticTail, |
| 18 | + StochasticTrapezoidalFins, |
| 19 | + StochasticParachute, |
| 20 | + StochasticRailButtons, |
| 21 | +) |
| 22 | +import datetime |
| 23 | + |
| 24 | +# If you are using Jupyter Notebooks, it is recommended to run the following line to make matplotlib plots which will be shown later interactive and higher quality. |
| 25 | +# |
| 26 | + |
| 27 | + |
| 28 | +# The `MonteCarlo` class allows to perform Monte Carlo simulations in a very simple way. |
| 29 | +# We just need to create an instance of the class, and then call the method `simulate()` to perform the simulations. |
| 30 | +# |
| 31 | +# The class has a lot of capabilities, we try to over as much as possible in this example. |
| 32 | +# We encourage you to check the documentation of the class to learn more about it. |
| 33 | +# |
| 34 | +# Additionally, you can check RocketPy's main reference for a better conceptual understanding |
| 35 | +# of the Monte Carlo Simulations: [RocketPy: Six Degree-of-Freedom Rocket Trajectory Simulator](<https://doi.org/10.1061/(ASCE)AS.1943-5525.0001331>) |
| 36 | + |
| 37 | +# ## Step 1: Creating the Inputs for the Simulations |
| 38 | +# |
| 39 | + |
| 40 | +# We will define rocketpy objects (e.g. Environment, SolidMotor, etc.) and use them to create the Monte Carlo simulation loop. |
| 41 | + |
| 42 | +# ### Environment |
| 43 | +# |
| 44 | + |
| 45 | +# Let's start by creating an Environment object, which will describe the atmospheric conditions for our launch site. |
| 46 | + |
| 47 | +env = Environment(latitude=39.389700, longitude=-8.288964, elevation=113) |
| 48 | + |
| 49 | +tomorrow = datetime.date.today() + datetime.timedelta(days=1) |
| 50 | + |
| 51 | +env.set_date((tomorrow.year, tomorrow.month, tomorrow.day, 12)) |
| 52 | + |
| 53 | +# As of the Atmospheric model, We will use an Ensemble model to a more accurate wind representation. |
| 54 | +# The Ensemble forecast is a different kind of forecast that uses multiple runs of the same model with slightly different initial conditions to represent the uncertainty in the forecast. |
| 55 | + |
| 56 | +env.set_atmospheric_model(type="Ensemble", file="GEFS") |
| 57 | + |
| 58 | +# For each rocketpy object, we will create a Stochastic object that works as an extension of the initial model, allow us to define what are the uncertainties of each input parameter. |
| 59 | + |
| 60 | +mc_env = StochasticEnvironment( |
| 61 | + environment=env, |
| 62 | + ensemble_member=list(range(env.num_ensemble_members)), |
| 63 | + wind_velocity_x_factor=(1.0, 0.33, "normal"), |
| 64 | + wind_velocity_y_factor=(1.0, 0.33, "normal"), |
| 65 | +) |
| 66 | + |
| 67 | +mc_env |
| 68 | + |
| 69 | +# In the example above, the |
| 70 | + |
| 71 | +# We can use the |
| 72 | + |
| 73 | +wind_x_at_1000m = [] |
| 74 | +for i in range(10): |
| 75 | + rnd_env = mc_env.create_object() |
| 76 | + wind_x_at_1000m.append(rnd_env.wind_velocity_x(1000)) |
| 77 | + |
| 78 | +# print(wind_x_at_1000m) |
| 79 | + |
| 80 | +# ### Motor |
| 81 | +# |
| 82 | + |
| 83 | +# Let's define the motor using the firs method. We will be using the data from the manufacturer, and following |
| 84 | +# the [RocketPy's documentation](https://docs.rocketpy.org/en/latest/user/index.html). |
| 85 | +# |
| 86 | + |
| 87 | +motor = SolidMotor( |
| 88 | + thrust_source="data/motors/Cesaroni_M1670.eng", |
| 89 | + dry_mass=1.815, |
| 90 | + dry_inertia=(0.125, 0.125, 0.002), |
| 91 | + nozzle_radius=33 / 1000, |
| 92 | + grain_number=5, |
| 93 | + grain_density=1815, |
| 94 | + grain_outer_radius=33 / 1000, |
| 95 | + grain_initial_inner_radius=15 / 1000, |
| 96 | + grain_initial_height=120 / 1000, |
| 97 | + grain_separation=5 / 1000, |
| 98 | + grains_center_of_mass_position=0.397, |
| 99 | + center_of_dry_mass_position=0.317, |
| 100 | + nozzle_position=0, |
| 101 | + burn_time=3.9, |
| 102 | + throat_radius=11 / 1000, |
| 103 | + coordinate_system_orientation="nozzle_to_combustion_chamber", |
| 104 | +) |
| 105 | +motor.total_impulse |
| 106 | + |
| 107 | +mc_motor = StochasticSolidMotor( |
| 108 | + solid_motor=motor, |
| 109 | + thrust_source=[ |
| 110 | + "data/motors/Cesaroni_M1670.eng", |
| 111 | + [[0, 6000], [1, 6000], [2, 6000], [3, 6000], [4, 6000]], |
| 112 | + Function([[0, 6000], [1, 6000], [2, 6000], [3, 6000], [4, 6000]]), |
| 113 | + ], |
| 114 | + burn_start_time=(0, 0.1), |
| 115 | + grains_center_of_mass_position=0.001, |
| 116 | + grain_density=50, |
| 117 | + grain_separation=1 / 1000, |
| 118 | + grain_initial_height=1 / 1000, |
| 119 | + grain_initial_inner_radius=0.375 / 1000, |
| 120 | + grain_outer_radius=0.375 / 1000, |
| 121 | + total_impulse=(6500, 1000), |
| 122 | + throat_radius=0.5 / 1000, |
| 123 | + nozzle_radius=0.5 / 1000, |
| 124 | + nozzle_position=0.001, |
| 125 | +) |
| 126 | +mc_motor |
| 127 | + |
| 128 | +total_impulse = [] |
| 129 | +for i in range(10): |
| 130 | + rnd_motor = mc_motor.create_object() |
| 131 | + total_impulse.append(rnd_motor.total_impulse) |
| 132 | + |
| 133 | +# print(total_impulse) |
| 134 | + |
| 135 | +# ### Rocket |
| 136 | +# |
| 137 | + |
| 138 | +rocket = Rocket( |
| 139 | + radius=127 / 2000, |
| 140 | + mass=14.426, |
| 141 | + inertia=(6.321, 6.321, 0.034), |
| 142 | + power_off_drag="data/calisto/powerOffDragCurve.csv", |
| 143 | + power_on_drag="data/calisto/powerOnDragCurve.csv", |
| 144 | + center_of_mass_without_motor=0, |
| 145 | + coordinate_system_orientation="tail_to_nose", |
| 146 | +) |
| 147 | + |
| 148 | +rail_buttons = rocket.set_rail_buttons( |
| 149 | + upper_button_position=0.0818, |
| 150 | + lower_button_position=-0.618, |
| 151 | + angular_position=45, |
| 152 | +) |
| 153 | + |
| 154 | +rocket.add_motor(motor, position=-1.255) |
| 155 | + |
| 156 | +nose_cone = rocket.add_nose(length=0.55829, kind="vonKarman", position=1.278) |
| 157 | + |
| 158 | +fin_set = rocket.add_trapezoidal_fins( |
| 159 | + n=4, |
| 160 | + root_chord=0.120, |
| 161 | + tip_chord=0.060, |
| 162 | + span=0.110, |
| 163 | + position=-1.04956, |
| 164 | + cant_angle=0.5, |
| 165 | + airfoil=("data/calisto/NACA0012-radians.csv", "radians"), |
| 166 | +) |
| 167 | + |
| 168 | +tail = rocket.add_tail( |
| 169 | + top_radius=0.0635, bottom_radius=0.0435, length=0.060, position=-1.194656 |
| 170 | +) |
| 171 | + |
| 172 | +# Additionally, we set parachutes for our Rocket, as well as the trigger functions for the deployment of such parachutes. |
| 173 | +# |
| 174 | + |
| 175 | +Main = rocket.add_parachute( |
| 176 | + "Main", |
| 177 | + cd_s=10.0, |
| 178 | + trigger=800, |
| 179 | + sampling_rate=105, |
| 180 | + lag=1.5, |
| 181 | + noise=(0, 8.3, 0.5), |
| 182 | +) |
| 183 | + |
| 184 | +Drogue = rocket.add_parachute( |
| 185 | + "Drogue", |
| 186 | + cd_s=1.0, |
| 187 | + trigger="apogee", |
| 188 | + sampling_rate=105, |
| 189 | + lag=1.5, |
| 190 | + noise=(0, 8.3, 0.5), |
| 191 | +) |
| 192 | + |
| 193 | +mc_rocket = StochasticRocket( |
| 194 | + rocket=rocket, |
| 195 | + radius=0.0127 / 2000, |
| 196 | + mass=(15.426, 0.5, "normal"), |
| 197 | + inertia_11=(6.321, 0), |
| 198 | + inertia_22=0.01, |
| 199 | + inertia_33=0.01, |
| 200 | + center_of_mass_without_motor=0, |
| 201 | +) |
| 202 | +mc_rocket |
| 203 | + |
| 204 | +mc_nose_cone = StochasticNoseCone( |
| 205 | + nosecone=nose_cone, |
| 206 | + length=0.001, |
| 207 | +) |
| 208 | + |
| 209 | +mc_fin_set = StochasticTrapezoidalFins( |
| 210 | + trapezoidal_fins=fin_set, |
| 211 | + root_chord=0.0005, |
| 212 | + tip_chord=0.0005, |
| 213 | + span=0.0005, |
| 214 | +) |
| 215 | + |
| 216 | +mc_tail = StochasticTail( |
| 217 | + tail=tail, |
| 218 | + top_radius=0.001, |
| 219 | + bottom_radius=0.001, |
| 220 | + length=0.001, |
| 221 | +) |
| 222 | + |
| 223 | +mc_rail_buttons = StochasticRailButtons( |
| 224 | + rail_buttons=rail_buttons, buttons_distance=0.001 |
| 225 | +) |
| 226 | + |
| 227 | +mc_main = StochasticParachute( |
| 228 | + parachute=Main, |
| 229 | + cd_s=0.1, |
| 230 | + lag=0.1, |
| 231 | +) |
| 232 | + |
| 233 | +mc_drogue = StochasticParachute( |
| 234 | + parachute=Drogue, |
| 235 | + cd_s=0.07, |
| 236 | + lag=0.2, |
| 237 | +) |
| 238 | + |
| 239 | +mc_rocket.add_motor(mc_motor, position=0.001) |
| 240 | +mc_rocket.add_nose(mc_nose_cone, position=(1.134, 0.001)) |
| 241 | +mc_rocket.add_trapezoidal_fins(mc_fin_set, position=(0.001, "normal")) |
| 242 | +mc_rocket.add_tail(mc_tail) |
| 243 | +mc_rocket.set_rail_buttons(mc_rail_buttons, lower_button_position=(0.001, "normal")) |
| 244 | +mc_rocket.add_parachute(mc_main) |
| 245 | +mc_rocket.add_parachute(mc_drogue) |
| 246 | + |
| 247 | +mc_rocket |
| 248 | + |
| 249 | +mc_rocket.last_rnd_dict |
| 250 | + |
| 251 | +# for i in range(10): |
| 252 | +# rnd_rocket = mc_rocket.create_object() |
| 253 | +# print(rnd_rocket.static_margin(0)) |
| 254 | + |
| 255 | +# |
| 256 | +# ### Flight |
| 257 | +# |
| 258 | + |
| 259 | +test_flight = Flight( |
| 260 | + rocket=rocket, |
| 261 | + environment=env, |
| 262 | + rail_length=5, |
| 263 | + inclination=84, |
| 264 | + heading=133, |
| 265 | +) |
| 266 | + |
| 267 | +mc_flight = StochasticFlight( |
| 268 | + flight=test_flight, |
| 269 | + inclination=(84.7, 1), |
| 270 | + heading=(53, 2), |
| 271 | +) |
| 272 | +mc_flight |
| 273 | + |
| 274 | +# And we can visualize the flight trajectory: |
| 275 | +# |
| 276 | + |
| 277 | +test_flight.plots.trajectory_3d() |
| 278 | + |
| 279 | +# ### Step 2: Starting the Monte Carlo Simulations |
| 280 | +# |
| 281 | + |
| 282 | +# First, let's invoke the `MonteCarlo` class, we only need a filename to initialize it. |
| 283 | +# The filename will be used either to save the results of the simulations or to load them |
| 284 | +# from a previous ran simulation. |
| 285 | +# |
| 286 | + |
| 287 | +test_dispersion = MonteCarlo( |
| 288 | + filename="monte_carlo_analysis_outputs/monte_carlo_class_example", |
| 289 | + environment=mc_env, |
| 290 | + rocket=mc_rocket, |
| 291 | + flight=mc_flight, |
| 292 | +) |
| 293 | + |
| 294 | +# TODO: add custom warning o when the rocket doesn't have a motors, parachute, or AeroSurface |
| 295 | + |
| 296 | +# Then, we can run the simulations using the method `MonteCarlo.simulate()`. |
| 297 | +# But before that, we need to set some simple parameters for the simulations. |
| 298 | +# We will set them by using a dictionary, which is one of the simplest way to do it. |
| 299 | +# |
| 300 | + |
| 301 | +# Finally, let's iterate over the simulations and export the data from each flight simulation! |
| 302 | +# |
| 303 | + |
| 304 | +test_dispersion.simulate(number_of_simulations=10, append=False, parallel=True) |
| 305 | + |
| 306 | +# ### Visualizing the results |
| 307 | +# |
| 308 | + |
| 309 | +# Now we finally have the results of our Monte Carlo simulations loaded! |
| 310 | +# Let's play with them. |
| 311 | +# |
| 312 | + |
| 313 | +# First, we can print numerical information regarding the results of the simulations. |
| 314 | +# |
| 315 | + |
| 316 | +# only need to import results if you did not run the simulations |
| 317 | +# test_dispersion.import_results() |
| 318 | + |
| 319 | +# test_dispersion.num_of_loaded_sims |
| 320 | + |
| 321 | +# test_dispersion.prints.all() |
| 322 | + |
| 323 | +# # Secondly, we can plot the results of the simulations. |
| 324 | +# # |
| 325 | + |
| 326 | +# test_dispersion.plots.ellipses(xlim=(-200, 3500), ylim=(-200, 3500)) |
| 327 | + |
| 328 | +# test_dispersion.plots.all() |
| 329 | + |
| 330 | +# Finally, one may also export to kml so it can be visualized in Google Earth |
| 331 | +# |
| 332 | + |
| 333 | +# test_dispersion.export_ellipses_to_kml( |
| 334 | +# filename="monte_carlo_analysis_outputs/monte_carlo_class_example.kml", |
| 335 | +# origin_lat=env.latitude, |
| 336 | +# origin_lon=env.longitude, |
| 337 | +# type="impact", |
| 338 | +# ) |
| 339 | + |
| 340 | + |
0 commit comments