Skip to content

Commit af03adf

Browse files
committed
Merge branch 'main' into symbolics
2 parents 0d0d449 + 0078822 commit af03adf

1 file changed

Lines changed: 348 additions & 4 deletions

File tree

docs/src/tutorial-continuation.md

Lines changed: 348 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,9 @@
11
# [Discrete continuation](@id tutorial-continuation)
22

3+
```@meta
4+
Draft = false
5+
```
6+
37
By using the warm start option, it is easy to implement a basic discrete continuation method, in which a sequence of problems is solved by using each solution as the initial guess for the next problem. This approach typically leads to faster and more reliable convergence than solving each problem with the same initial guess and is particularly useful for problems that require a good initial guess to converge.
48

59
## Continuation on parametric OCP
@@ -128,20 +132,18 @@ function goddard_problem(Tmax)
128132
129133
tf ∈ R, variable
130134
t ∈ [0, tf], time
131-
x ∈ R^3, state
135+
x = (r, v, m) ∈ R^3, state
132136
u ∈ R, control
133137
134138
0.01 ≤ tf ≤ Inf
135139
136-
r = x[1]
137-
v = x[2]
138-
m = x[3]
139140
x(0) == x0
140141
m(tf) == mf
141142
r0 ≤ r(t) ≤ r0 + 0.1
142143
v0 ≤ v(t) ≤ vmax
143144
mf ≤ m(t) ≤ m0
144145
0 ≤ u(t) ≤ 1
146+
145147
ẋ(t) == F0(x(t)) + u(t) * F1(x(t), Tmax)
146148
147149
r(tf) → max
@@ -212,6 +214,348 @@ end
212214
plot(plt_sol; size=(800, 800), leftmargin=5mm)
213215
```
214216

217+
### Animation
218+
219+
The animation below shows three rockets launching simultaneously with different maximum thrust values (Tmax = 3, 2, 1). The rockets start at the same time (t=0) but reach their final altitudes at different times according to their optimal trajectories. Each rocket displays:
220+
221+
- **Altitude** (above the rocket)
222+
- **Velocity gauge** (horizontal bar, green→red)
223+
- **Fuel gauge** (vertical bar, showing remaining fuel m(t) - mf)
224+
- **Flame** (visible when thrust u(t) > 0)
225+
226+
```@setup main-cont
227+
# Extract data for the 3 selected solutions (Tmax ≈ 3, 2, 1)
228+
selected_sols = [sols[idx] for idx in indices]
229+
selected_Tmax = [round(data.Tmax[idx], digits=2) for idx in indices]
230+
231+
# Extract time grids and state/control data for each solution
232+
solutions_data = []
233+
for (sol, Tmax_val) in zip(selected_sols, selected_Tmax)
234+
t_grid = time_grid(sol)
235+
X = state(sol).(t_grid)
236+
u = control(sol).(t_grid)
237+
238+
# Extract r, v, m from state
239+
r = [x[1] for x in X]
240+
v = [x[2] for x in X]
241+
m = [x[3] for x in X]
242+
243+
push!(solutions_data, (t=t_grid, r=r, v=v, m=m, u=u))
244+
end
245+
246+
# Serialize to JSON for JavaScript injection
247+
json_t1 = "[" * join(string.(solutions_data[1].t), ",") * "]"
248+
json_r1 = "[" * join(string.(solutions_data[1].r), ",") * "]"
249+
json_v1 = "[" * join(string.(solutions_data[1].v), ",") * "]"
250+
json_m1 = "[" * join(string.(solutions_data[1].m), ",") * "]"
251+
json_u1 = "[" * join(string.(solutions_data[1].u), ",") * "]"
252+
253+
json_t2 = "[" * join(string.(solutions_data[2].t), ",") * "]"
254+
json_r2 = "[" * join(string.(solutions_data[2].r), ",") * "]"
255+
json_v2 = "[" * join(string.(solutions_data[2].v), ",") * "]"
256+
json_m2 = "[" * join(string.(solutions_data[2].m), ",") * "]"
257+
json_u2 = "[" * join(string.(solutions_data[2].u), ",") * "]"
258+
259+
json_t3 = "[" * join(string.(solutions_data[3].t), ",") * "]"
260+
json_r3 = "[" * join(string.(solutions_data[3].r), ",") * "]"
261+
json_v3 = "[" * join(string.(solutions_data[3].v), ",") * "]"
262+
json_m3 = "[" * join(string.(solutions_data[3].m), ",") * "]"
263+
json_u3 = "[" * join(string.(solutions_data[3].u), ",") * "]"
264+
265+
# Inject real Tmax values for labels
266+
json_Tmax1 = string(selected_Tmax[1])
267+
json_Tmax2 = string(selected_Tmax[2])
268+
json_Tmax3 = string(selected_Tmax[3])
269+
270+
# Inject m0 and mf values
271+
json_m0 = string(m0)
272+
json_mf = string(mf)
273+
274+
# Inject r0 value
275+
json_r0 = string(r0)
276+
277+
# Calculate global max time for animation loop
278+
t_max_global = max(solutions_data[1].t[end], solutions_data[2].t[end], solutions_data[3].t[end])
279+
280+
# Define RawHTML wrapper for Documenter
281+
struct RawHTML
282+
raw::String
283+
end
284+
Base.show(io::IO, ::MIME"text/html", h::RawHTML) = print(io, h.raw)
285+
286+
html_anim = """
287+
<div style="display: flex; justify-content: center; margin: 20px 0;">
288+
<canvas id="goddardCanvas" width="900" height="650"
289+
style="border:1px solid #ddd; border-radius: 8px; max-width: 100%;">
290+
</canvas>
291+
</div>
292+
293+
<script>
294+
(function() {
295+
// Data injected from Julia
296+
const t1 = $json_t1, r1 = $json_r1, v1 = $json_v1, m1 = $json_m1, u1 = $json_u1;
297+
const t2 = $json_t2, r2 = $json_r2, v2 = $json_v2, m2 = $json_m2, u2 = $json_u2;
298+
const t3 = $json_t3, r3 = $json_r3, v3 = $json_v3, m3 = $json_m3, u3 = $json_u3;
299+
300+
const Tmax1 = $json_Tmax1, Tmax2 = $json_Tmax2, Tmax3 = $json_Tmax3;
301+
302+
const t_max_global = $t_max_global;
303+
const m0 = $json_m0, mf = $json_mf;
304+
const r0 = $json_r0;
305+
306+
const canvas = document.getElementById('goddardCanvas');
307+
const ctx = canvas.getContext('2d');
308+
309+
// Detect dark/light mode
310+
const isDark = window.matchMedia && window.matchMedia('(prefers-color-scheme: dark)').matches;
311+
312+
// Color palettes
313+
const colors = isDark ? {
314+
canvas_bg: '#1e1e2e',
315+
ground: '#4a4a4a',
316+
gauge_bg: '#3a3a3a',
317+
gauge_high: '#2ecc71',
318+
gauge_mid: '#f39c12',
319+
gauge_low: '#e74c3c',
320+
text: '#ecf0f1',
321+
text_secondary: '#bdc3c7',
322+
flame_outer: '#e67e22',
323+
flame_inner: '#f1c40f',
324+
window_fill: 'rgba(255,255,255,0.3)',
325+
window_stroke: 'rgba(255,255,255,0.6)',
326+
progress: '#4063D8'
327+
} : {
328+
canvas_bg: '#fafafa',
329+
ground: '#ccc',
330+
gauge_bg: '#ecf0f1',
331+
gauge_high: '#2ecc71',
332+
gauge_mid: '#f39c12',
333+
gauge_low: '#e74c3c',
334+
text: '#2c3e50',
335+
text_secondary: '#7f8c8d',
336+
flame_outer: '#e67e22',
337+
flame_inner: '#f1c40f',
338+
window_fill: 'rgba(255,255,255,0.55)',
339+
window_stroke: 'rgba(255,255,255,0.85)',
340+
progress: '#4063D8'
341+
};
342+
343+
// Animation duration in seconds (independent of real problem time)
344+
const animation_duration = 10.0;
345+
346+
const zone_width = canvas.width / 3;
347+
const top_margin = 50; // Space for time label at top
348+
const bottom_margin = 120; // Space for ground line + v gauge + Tmax label
349+
const baseline = canvas.height - bottom_margin;
350+
351+
// Scale for altitude (r goes from r0 to ~r0+0.1)
352+
const r_min = r0;
353+
const r_max_all = Math.max(...r1, ...r2, ...r3);
354+
const r_range = r_max_all - r_min;
355+
const altitude_scale = (baseline - top_margin) / r_range;
356+
357+
// Scale for velocity (v goes from 0 to ~0.1)
358+
const v_max_all = Math.max(...v1, ...v2, ...v3);
359+
360+
let start_time = null;
361+
362+
function interpolate(t, data_t, data_y) {
363+
// Find interval
364+
let i = 0;
365+
while(i < data_t.length - 1 && data_t[i + 1] < t) { i++; }
366+
367+
if (i >= data_t.length - 1) {
368+
return data_y[data_y.length - 1]; // Return final value if past end
369+
}
370+
371+
const dt = data_t[i+1] - data_t[i];
372+
const alpha = (dt > 0) ? (t - data_t[i]) / dt : 0;
373+
return data_y[i] + alpha * (data_y[i+1] - data_y[i]);
374+
}
375+
376+
function drawRocket(zone_idx, t, r, v, m, u, color, Tmax) {
377+
const zone_x = zone_idx * zone_width;
378+
const center_x = zone_x + zone_width / 2;
379+
380+
// Rocket position
381+
const rocket_y = baseline - altitude_scale * (r - r_min);
382+
383+
// Draw ground line
384+
ctx.beginPath();
385+
ctx.moveTo(zone_x + 20, baseline + 20);
386+
ctx.lineTo(zone_x + zone_width - 20, baseline + 20);
387+
ctx.lineWidth = 2;
388+
ctx.strokeStyle = colors.ground;
389+
ctx.stroke();
390+
391+
// Draw fuel gauge (vertical bar to the right of rocket)
392+
const fuel_x = center_x + 80;
393+
const fuel_height = 200;
394+
const fuel_y = baseline - fuel_height;
395+
const fuel_level = (m - mf) / (m0 - mf); // 0 to 1
396+
397+
// Fuel gauge background
398+
ctx.fillStyle = colors.gauge_bg;
399+
ctx.fillRect(fuel_x - 10, fuel_y, 20, fuel_height);
400+
401+
// Fuel gauge fill
402+
const fuel_fill_height = fuel_level * fuel_height;
403+
ctx.fillStyle = fuel_level > 0.5 ? colors.gauge_high : (fuel_level > 0.2 ? colors.gauge_mid : colors.gauge_low);
404+
ctx.fillRect(fuel_x - 10, fuel_y + (fuel_height - fuel_fill_height), 20, fuel_fill_height);
405+
406+
// Fuel gauge label
407+
ctx.fillStyle = colors.text_secondary;
408+
ctx.font = '14px Arial';
409+
ctx.textAlign = 'center';
410+
ctx.fillText('Fuel', fuel_x, fuel_y - 8);
411+
412+
// Draw rocket: nose cone + body + fins + window
413+
const rocket_bottom = rocket_y + 18;
414+
ctx.fillStyle = color;
415+
416+
// Nose cone
417+
ctx.beginPath();
418+
ctx.moveTo(center_x, rocket_y - 38);
419+
ctx.lineTo(center_x - 12, rocket_y - 10);
420+
ctx.lineTo(center_x + 12, rocket_y - 10);
421+
ctx.closePath();
422+
ctx.fill();
423+
424+
// Body
425+
ctx.fillRect(center_x - 12, rocket_y - 10, 24, 28);
426+
427+
// Left fin
428+
ctx.beginPath();
429+
ctx.moveTo(center_x - 12, rocket_y + 6);
430+
ctx.lineTo(center_x - 24, rocket_y + 22);
431+
ctx.lineTo(center_x - 12, rocket_y + 18);
432+
ctx.closePath();
433+
ctx.fill();
434+
435+
// Right fin
436+
ctx.beginPath();
437+
ctx.moveTo(center_x + 12, rocket_y + 6);
438+
ctx.lineTo(center_x + 24, rocket_y + 22);
439+
ctx.lineTo(center_x + 12, rocket_y + 18);
440+
ctx.closePath();
441+
ctx.fill();
442+
443+
// Window
444+
ctx.beginPath();
445+
ctx.arc(center_x, rocket_y + 2, 5, 0, 2 * Math.PI);
446+
ctx.fillStyle = colors.window_fill;
447+
ctx.fill();
448+
ctx.strokeStyle = colors.window_stroke;
449+
ctx.lineWidth = 1.5;
450+
ctx.stroke();
451+
452+
// Draw flame if thrust > 0.1
453+
if (u > 0.1) {
454+
const flame_size = 10 + 20 * u;
455+
ctx.fillStyle = colors.flame_outer;
456+
ctx.beginPath();
457+
ctx.moveTo(center_x - 10, rocket_bottom);
458+
ctx.lineTo(center_x, rocket_bottom + flame_size);
459+
ctx.lineTo(center_x + 10, rocket_bottom);
460+
ctx.closePath();
461+
ctx.fill();
462+
463+
// Inner flame
464+
ctx.fillStyle = colors.flame_inner;
465+
ctx.beginPath();
466+
ctx.moveTo(center_x - 5, rocket_bottom);
467+
ctx.lineTo(center_x, rocket_bottom + flame_size * 0.6);
468+
ctx.lineTo(center_x + 5, rocket_bottom);
469+
ctx.closePath();
470+
ctx.fill();
471+
}
472+
473+
// Draw altitude text
474+
ctx.fillStyle = colors.text;
475+
ctx.font = 'bold 14px Arial';
476+
ctx.textAlign = 'center';
477+
ctx.fillText('r = ' + r.toFixed(4), center_x, rocket_y - 40);
478+
479+
// Draw velocity gauge (horizontal bar)
480+
const v_gauge_width = 60;
481+
const v_gauge_height = 8;
482+
const v_gauge_x = center_x - v_gauge_width / 2;
483+
const v_gauge_y = rocket_y + 50;
484+
const v_level = Math.abs(v) / v_max_all;
485+
486+
// Velocity gauge background
487+
ctx.fillStyle = colors.gauge_bg;
488+
ctx.fillRect(v_gauge_x, v_gauge_y, v_gauge_width, v_gauge_height);
489+
490+
// Velocity gauge fill
491+
ctx.fillStyle = v_level > 0.7 ? colors.gauge_low : (v_level > 0.4 ? colors.gauge_mid : colors.gauge_high);
492+
ctx.fillRect(v_gauge_x, v_gauge_y, v_level * v_gauge_width, v_gauge_height);
493+
494+
// Velocity label
495+
ctx.fillStyle = colors.text_secondary;
496+
ctx.font = '12px Arial';
497+
ctx.fillText('v', center_x, v_gauge_y + 20);
498+
499+
// Draw Tmax label at bottom
500+
ctx.fillStyle = colors.text;
501+
ctx.font = 'bold 16px Arial';
502+
ctx.fillText('Tmax = ' + Tmax, center_x, baseline + 90);
503+
}
504+
505+
function draw(time) {
506+
if (!start_time) start_time = time;
507+
508+
const elapsed = (time - start_time) / 1000.0;
509+
const progress = (elapsed % animation_duration) / animation_duration;
510+
const t_anim = progress * t_max_global;
511+
512+
ctx.fillStyle = colors.canvas_bg;
513+
ctx.fillRect(0, 0, canvas.width, canvas.height);
514+
515+
// Draw time indicator at top
516+
ctx.fillStyle = colors.text;
517+
ctx.font = 'bold 18px Arial';
518+
ctx.textAlign = 'center';
519+
ctx.fillText('t = ' + t_anim.toFixed(3) + ' s', canvas.width / 2, 30);
520+
521+
// Interpolate and draw each rocket
522+
const r1_cur = interpolate(t_anim, t1, r1);
523+
const v1_cur = interpolate(t_anim, t1, v1);
524+
const m1_cur = interpolate(t_anim, t1, m1);
525+
const u1_cur = interpolate(t_anim, t1, u1);
526+
drawRocket(0, t_anim, r1_cur, v1_cur, m1_cur, u1_cur, '#CB3C33', Tmax1);
527+
528+
const r2_cur = interpolate(t_anim, t2, r2);
529+
const v2_cur = interpolate(t_anim, t2, v2);
530+
const m2_cur = interpolate(t_anim, t2, m2);
531+
const u2_cur = interpolate(t_anim, t2, u2);
532+
drawRocket(1, t_anim, r2_cur, v2_cur, m2_cur, u2_cur, '#389826', Tmax2);
533+
534+
const r3_cur = interpolate(t_anim, t3, r3);
535+
const v3_cur = interpolate(t_anim, t3, v3);
536+
const m3_cur = interpolate(t_anim, t3, m3);
537+
const u3_cur = interpolate(t_anim, t3, u3);
538+
drawRocket(2, t_anim, r3_cur, v3_cur, m3_cur, u3_cur, '#9558B2', Tmax3);
539+
540+
// Draw progress bar at bottom
541+
const bar_height = 4;
542+
ctx.fillStyle = colors.gauge_bg;
543+
ctx.fillRect(0, canvas.height - bar_height, canvas.width, bar_height);
544+
ctx.fillStyle = colors.progress;
545+
ctx.fillRect(0, canvas.height - bar_height, canvas.width * progress, bar_height);
546+
547+
requestAnimationFrame(draw);
548+
}
549+
requestAnimationFrame(draw);
550+
})();
551+
</script>
552+
"""
553+
```
554+
555+
```@example main-cont
556+
RawHTML(html_anim) # hide
557+
```
558+
215559
## Best practices and limitations
216560

217561
The examples above illustrate a basic discrete continuation method. Here are some important considerations:

0 commit comments

Comments
 (0)