UNPKG

mathjs/examples/browser/rocket_trajectory_optimization.html

Version:

1<!DOCTYPE html>
2<html lang="en">
3<head>
<meta charset="utf-8">
<title>math.js | rocket trajectory optimization</title>
6
<script src="../../dist/math.js"></script>
<script src="https://cdnjs.cloudflare.com/ajax/libs/Chart.js/2.5.0/Chart.min.js"></script>
9
<style>
  body {
    font-family: sans-serif;
  }
</style>
15</head>
16<body>
<h1>Rocket trajectory optimization</h1>
<p>
  This example simulates the ascent stage of the Apollo Lunar Module modeled using a system of ordinary differential equations.
</p>
21
<canvas id=canvas1 width=1600 height=400></canvas>
<canvas id=canvas2 width=800 height=400></canvas>
24
<script>
26
  function ndsolve(f, x0, dt, tmax) {
    const n = f.size()[0]  // Number of variables
    const x = x0.clone()   // Current values of variables
    const dxdt = []        // Temporary variable to hold time-derivatives
    const result = []      // Contains entire solution
32
    const nsteps = math.divide(tmax, dt)   // Number of time steps
    for(let i=0; i<nsteps; i++) {
      // Compute derivatives
      for(let j=0; j<n; j++) {
        dxdt[j] = f.get([j]).apply(null, x.toArray())
      }
      // Euler method to compute next time step
      for(let j=0; j<n; j++) {
        x.set([j], math.add(x.get([j]), math.multiply(dxdt[j], dt)))
      }
      result.push(x.clone())
    }
45
    return math.matrix(result)
  }
48
  // Import the numerical ODE solver
  math.import({ndsolve:ndsolve})
51
  // Create a math.js context for our simulation. Everything else occurs in the context of the expression parser!
53
  const sim = math.parser()
55
  sim.evaluate("G = 6.67408e-11 m^3 kg^-1 s^-2")  // Gravitational constant
  sim.evaluate("mbody = 5.972e24 kg")             // Mass of Earth
  sim.evaluate("mu = G * mbody")
  sim.evaluate("dt = 1.0 s")                      // Simulation timestep
  sim.evaluate("tfinal = 162 s")                  // Simulation duration
  sim.evaluate("T = 1710000 lbf * 0.9")           // Engine thrust
  sim.evaluate("g0 = 9.80665 m/s^2")              // Standard gravity: used for calculating prop consumption (dmdt)
  sim.evaluate("isp = 290 s")                     // Specific impulse
  sim.evaluate("gamma0 = 89.99883 deg")           // Initial pitch angle (90 deg is vertical)
  sim.evaluate("r0 = 6378.1370 km")               // Equatorial radius of Earth
  sim.evaluate("v0 = 10 m/s")                     // Initial velocity (must be non-zero because ODE is ill-conditioned)
  sim.evaluate("phi0 = 0 deg")                    // Initial orbital reference angle
  sim.evaluate("m0 = 1207920 lbm + 30000 lbm")    // Initial mass of rocket and fuel
69
  // Define the equations of motion. It is important to maintain the same argument order for each of these functions.
  sim.evaluate("drdt(r, v, m, phi, gamma) = v sin(gamma)")
  sim.evaluate("dvdt(r, v, m, phi, gamma) = -mu / r^2 * sin(gamma) + T / m")
  sim.evaluate("dmdt(r, v, m, phi, gamma) = -T/g0/isp")
  sim.evaluate("dphidt(r, v, m, phi, gamma) = v/r * cos(gamma) * rad")
  sim.evaluate("dgammadt(r, v, m, phi, gamma) = (1/r * (v - mu / (r v)) * cos(gamma)) * rad")
76
  // Again, remember to maintain the same variable order in the call to ndsolve.
  sim.evaluate("result_stage1 = ndsolve([drdt, dvdt, dmdt, dphidt, dgammadt], [r0, v0, m0, phi0, gamma0], dt, tfinal)")
79
  // Reset initial conditions for interstage flight
  sim.evaluate("T = 0 lbf")
  sim.evaluate("tfinal = 12 s")
  sim.evaluate("x = flatten(result_stage1[result_stage1.size()[1],:])")
  sim.evaluate("result_interstage = ndsolve([drdt, dvdt, dmdt, dphidt, dgammadt], x, dt, tfinal)")
85
  console.log(sim.evaluate("result_interstage[result_interstage.size()[1],3]").toString())
87
  // Reset initial conditions for stage 2 flight
  sim.evaluate("T = 210000 lbf")
  sim.evaluate("isp = 348 s")
  sim.evaluate("tfinal = 397 s")
  sim.evaluate("x = flatten(result_interstage[result_interstage.size()[1],:])")
  sim.evaluate("x[3] = 273600 lbm")  // Lighten the rocket a bit since we discarded the first stage
  sim.evaluate("result_stage2 = ndsolve([drdt, dvdt, dmdt, dphidt, dgammadt], x, dt, tfinal)")
95
  // Reset initial conditions for unpowered flight
  sim.evaluate("T = 0 lbf")
  sim.evaluate("tfinal = 60 s")
  sim.evaluate("x = flatten(result_stage2[result_stage2.size()[1],:])")
  sim.evaluate("result_unpowered = ndsolve([drdt, dvdt, dmdt, dphidt, dgammadt], x, dt, tfinal)")
101
102
103
  // Extract the useful information from the results so it can be plotted
  const data_stage1 =     sim.evaluate("transpose(concat( transpose(    result_stage1[:,4] - phi0) * r0 / rad / km, (    transpose(result_stage1[:,1]) - r0) / km, 1 ))").toArray().map(function(e) { return {x: e[0], y: e[1]} })
  const data_interstage = sim.evaluate("transpose(concat( transpose(result_interstage[:,4] - phi0) * r0 / rad / km, (transpose(result_interstage[:,1]) - r0) / km, 1 ))").toArray().map(function(e) { return {x: e[0], y: e[1]} })
  const data_stage2 =     sim.evaluate("transpose(concat( transpose(    result_stage2[:,4] - phi0) * r0 / rad / km, (    transpose(result_stage2[:,1]) - r0) / km, 1 ))").toArray().map(function(e) { return {x: e[0], y: e[1]} })
  const data_unpowered =  sim.evaluate("transpose(concat( transpose( result_unpowered[:,4] - phi0) * r0 / rad / km, ( transpose(result_unpowered[:,1]) - r0) / km, 1 ))").toArray().map(function(e) { return {x: e[0], y: e[1]} })
109
  window['chart'] = new Chart(document.getElementById('canvas1'), {
    type: 'line',
    data: {
      datasets: [{
        label: "Stage 1",
        data: data_stage1,
        fill: false,
        borderColor: "red",
        pointRadius: 0
      }, {
        label: "Interstage",
        data: data_interstage,
        fill: false,
        borderColor: "green",
        pointRadius: 0
      }, {
        label: "Stage 2",
        data: data_stage2,
        fill: false,
        borderColor: "orange",
        pointRadius: 0
      }, {
        label: "Unpowered",
        data: data_unpowered,
        fill: false,
        borderColor: "blue",
        pointRadius: 0
      }]
    },
    options: {
      scales: {
        xAxes: [{
          type: 'linear',
          position: 'bottom'
        }]
      }
    }
147
  })
149
</script>
151</body>
152</html>

1	`<!DOCTYPE html>`
2	`<html lang="en">`
3	`<head>`
4	`<meta charset="utf-8">`
5	`<title>math.js \| rocket trajectory optimization</title>`
6
7	`<script src="../../dist/math.js"></script>`
8	`<script src="https://cdnjs.cloudflare.com/ajax/libs/Chart.js/2.5.0/Chart.min.js"></script>`
9
10	`<style>`
11	`body {`
12	`font-family: sans-serif;`
13	`}`
14	`</style>`
15	`</head>`
16	`<body>`
17	`<h1>Rocket trajectory optimization</h1>`
18	`<p>`
19	`This example simulates the ascent stage of the Apollo Lunar Module modeled using a system of ordinary differential equations.`
20	`</p>`
21
22	`<canvas id=canvas1 width=1600 height=400></canvas>`
23	`<canvas id=canvas2 width=800 height=400></canvas>`
24
25	`<script>`
26
27	`function ndsolve(f, x0, dt, tmax) {`
28	`const n = f.size()[0] // Number of variables`
29	`const x = x0.clone() // Current values of variables`
30	`const dxdt = [] // Temporary variable to hold time-derivatives`
31	`const result = [] // Contains entire solution`
32
33	`const nsteps = math.divide(tmax, dt) // Number of time steps`
34	`for(let i=0; i<nsteps; i++) {`
35	`// Compute derivatives`
36	`for(let j=0; j<n; j++) {`
37	`dxdt[j] = f.get([j]).apply(null, x.toArray())`
38	`}`
39	`// Euler method to compute next time step`
40	`for(let j=0; j<n; j++) {`
41	`x.set([j], math.add(x.get([j]), math.multiply(dxdt[j], dt)))`
42	`}`
43	`result.push(x.clone())`
44	`}`
45
46	`return math.matrix(result)`
47	`}`
48
49	`// Import the numerical ODE solver`
50	`math.import({ndsolve:ndsolve})`
51
52	`// Create a math.js context for our simulation. Everything else occurs in the context of the expression parser!`
53
54	`const sim = math.parser()`
55
56	`sim.evaluate("G = 6.67408e-11 m^3 kg^-1 s^-2") // Gravitational constant`
57	`sim.evaluate("mbody = 5.972e24 kg") // Mass of Earth`
58	`sim.evaluate("mu = G * mbody")`
59	`sim.evaluate("dt = 1.0 s") // Simulation timestep`
60	`sim.evaluate("tfinal = 162 s") // Simulation duration`
61	`sim.evaluate("T = 1710000 lbf * 0.9") // Engine thrust`
62	`sim.evaluate("g0 = 9.80665 m/s^2") // Standard gravity: used for calculating prop consumption (dmdt)`
63	`sim.evaluate("isp = 290 s") // Specific impulse`
64	`sim.evaluate("gamma0 = 89.99883 deg") // Initial pitch angle (90 deg is vertical)`
65	`sim.evaluate("r0 = 6378.1370 km") // Equatorial radius of Earth`
66	`sim.evaluate("v0 = 10 m/s") // Initial velocity (must be non-zero because ODE is ill-conditioned)`
67	`sim.evaluate("phi0 = 0 deg") // Initial orbital reference angle`
68	`sim.evaluate("m0 = 1207920 lbm + 30000 lbm") // Initial mass of rocket and fuel`
69
70	`// Define the equations of motion. It is important to maintain the same argument order for each of these functions.`
71	`sim.evaluate("drdt(r, v, m, phi, gamma) = v sin(gamma)")`
72	`sim.evaluate("dvdt(r, v, m, phi, gamma) = -mu / r^2 * sin(gamma) + T / m")`
73	`sim.evaluate("dmdt(r, v, m, phi, gamma) = -T/g0/isp")`
74	`sim.evaluate("dphidt(r, v, m, phi, gamma) = v/r * cos(gamma) * rad")`
75	`sim.evaluate("dgammadt(r, v, m, phi, gamma) = (1/r * (v - mu / (r v)) * cos(gamma)) * rad")`
76
77	`// Again, remember to maintain the same variable order in the call to ndsolve.`
78	`sim.evaluate("result_stage1 = ndsolve([drdt, dvdt, dmdt, dphidt, dgammadt], [r0, v0, m0, phi0, gamma0], dt, tfinal)")`
79
80	`// Reset initial conditions for interstage flight`
81	`sim.evaluate("T = 0 lbf")`
82	`sim.evaluate("tfinal = 12 s")`
83	`sim.evaluate("x = flatten(result_stage1[result_stage1.size()[1],:])")`
84	`sim.evaluate("result_interstage = ndsolve([drdt, dvdt, dmdt, dphidt, dgammadt], x, dt, tfinal)")`
85
86	`console.log(sim.evaluate("result_interstage[result_interstage.size()[1],3]").toString())`
87
88	`// Reset initial conditions for stage 2 flight`
89	`sim.evaluate("T = 210000 lbf")`
90	`sim.evaluate("isp = 348 s")`
91	`sim.evaluate("tfinal = 397 s")`
92	`sim.evaluate("x = flatten(result_interstage[result_interstage.size()[1],:])")`
93	`sim.evaluate("x[3] = 273600 lbm") // Lighten the rocket a bit since we discarded the first stage`
94	`sim.evaluate("result_stage2 = ndsolve([drdt, dvdt, dmdt, dphidt, dgammadt], x, dt, tfinal)")`
95
96	`// Reset initial conditions for unpowered flight`
97	`sim.evaluate("T = 0 lbf")`
98	`sim.evaluate("tfinal = 60 s")`
99	`sim.evaluate("x = flatten(result_stage2[result_stage2.size()[1],:])")`
100	`sim.evaluate("result_unpowered = ndsolve([drdt, dvdt, dmdt, dphidt, dgammadt], x, dt, tfinal)")`
101
102
103
104	`// Extract the useful information from the results so it can be plotted`
105	`const data_stage1 = sim.evaluate("transpose(concat( transpose( result_stage1[:,4] - phi0) * r0 / rad / km, ( transpose(result_stage1[:,1]) - r0) / km, 1 ))").toArray().map(function(e) { return {x: e[0], y: e[1]} })`
106	`const data_interstage = sim.evaluate("transpose(concat( transpose(result_interstage[:,4] - phi0) * r0 / rad / km, (transpose(result_interstage[:,1]) - r0) / km, 1 ))").toArray().map(function(e) { return {x: e[0], y: e[1]} })`
107	`const data_stage2 = sim.evaluate("transpose(concat( transpose( result_stage2[:,4] - phi0) * r0 / rad / km, ( transpose(result_stage2[:,1]) - r0) / km, 1 ))").toArray().map(function(e) { return {x: e[0], y: e[1]} })`
108	`const data_unpowered = sim.evaluate("transpose(concat( transpose( result_unpowered[:,4] - phi0) * r0 / rad / km, ( transpose(result_unpowered[:,1]) - r0) / km, 1 ))").toArray().map(function(e) { return {x: e[0], y: e[1]} })`
109
110	`window['chart'] = new Chart(document.getElementById('canvas1'), {`
111	`type: 'line',`
112	`data: {`
113	`datasets: [{`
114	`label: "Stage 1",`
115	`data: data_stage1,`
116	`fill: false,`
117	`borderColor: "red",`
118	`pointRadius: 0`
119	`}, {`
120	`label: "Interstage",`
121	`data: data_interstage,`
122	`fill: false,`
123	`borderColor: "green",`
124	`pointRadius: 0`
125	`}, {`
126	`label: "Stage 2",`
127	`data: data_stage2,`
128	`fill: false,`
129	`borderColor: "orange",`
130	`pointRadius: 0`
131	`}, {`
132	`label: "Unpowered",`
133	`data: data_unpowered,`
134	`fill: false,`
135	`borderColor: "blue",`
136	`pointRadius: 0`
137	`}]`
138	`},`
139	`options: {`
140	`scales: {`
141	`xAxes: [{`
142	`type: 'linear',`
143	`position: 'bottom'`
144	`}]`
145	`}`
146	`}`
147
148	`})`
149
150	`</script>`
151	`</body>`
152	`</html>`