가장 멀리 나는 속도: 1차원 최적화

둘째 주 월요일 아침, 민재 팀장이 비행 기록 화면을 띄운 채 관제실로 들어옵니다. "섬 보건소까지 21 km인데, 19 km쯤에서 배터리 경고가 떠서 돌려보냈어요. 배터리를 더 큰 걸로 바꿔야 할까요?"

서윤이 먼저 답합니다. "빨리 날면 빨리 도착하니까, 속도를 더 올리면 배터리를 덜 쓰지 않을까요?" 그러다 고개를 갸웃합니다. "아니다, 천천히 날면 모터가 힘을 덜 쓰니까 오히려 느리게 가야 하나?"

오 박사가 화이트보드에 언덕 모양 곡선을 하나 그립니다. "둘 다 반만 맞아. 너무 느리면 몸을 띄우는 데 힘을 다 쓰고, 너무 빠르면 공기를 밀어내는 데 힘을 다 쓰지. 그 사이 어딘가에 꼭대기가 있어. 꼭대기를 찾는 것도 한 걸음씩 가면 돼."

이번 교시의 주제는 최적화(optimization), 그중에서도 변수가 하나뿐인 1차원 최적화입니다. 4·5교시에서 "f(x) = 0인 곳"을 찾았다면, 이번에는 "f(x)가 가장 큰 곳"을 찾습니다. 두 문제는 생각보다 가깝고, 또 생각보다 다릅니다.

6-1. 속도마다 거리 계산하기: 모델과 격자 탐색 (8분)

드론이 속도 $v$로 수평 비행할 때 필요한 전력 $P(v)$(와트, W)는 크게 세 덩어리로 나눌 수 있습니다.

첫째 항은 프로펠러가 공기를 아래로 밀어 몸을 띄우는 데 드는 전력, 유도 전력(induced power) 입니다. 빨리 날수록 1초에 지나가는 공기가 많아서, 같은 무게를 띄우는 데 힘이 덜 듭니다. 그래서 $v$에 반비례합니다. 둘째 항은 공기 저항을 이기는 전력입니다. 1교시의 항력이 $v^2$에 비례했고, 전력은 힘 × 속도이므로 $v^3$에 비례합니다. 셋째 항은 비행 제어기, 카메라, 통신 장비가 늘 쓰는 전력입니다.

마지막 가정이 중요합니다. 이 식은 $v \to 0$이면 전력이 무한대가 되는데, 실제 드론은 제자리에 떠 있을 수 있습니다. 모델은 쓰려는 범위 안에서만 믿습니다.

배터리 에너지를 전력으로 나누면 비행 시간이고, 비행 시간에 속도를 곱하면 거리입니다. 단위를 맞추면 이렇습니다. Wh ÷ W = 시간(h)이므로 분으로는 60을 곱하고, m/s × h를 km로 바꾸려면 3600/1000 = 3.6을 곱합니다.

먼저 2 m/s 간격으로 표를 만들어 봅시다.

>>> E, A, B, C = 400, 4500, 0.14, 150
>>> def power(v):
...     return A / v + B * v**3 + C
>>> def flight_min(v):
...     return 60 * E / power(v)
>>> def range_km(v):
...     return 3.6 * v * E / power(v)
>>> v = np.arange(8, 23, 2)
>>> v
array([ 8, 10, 12, 14, 16, 18, 20, 22])
>>> np.round(power(v))
array([ 784.,  740.,  767.,  856., 1005., 1216., 1495., 1845.])
>>> np.round(flight_min(v), 1)
array([30.6, 32.4, 31.3, 28.1, 23.9, 19.7, 16.1, 13. ])
>>> np.round(range_km(v), 1)
array([14.7, 19.5, 22.5, 23.6, 22.9, 21.3, 19.3, 17.2])

20 m/s에서는 19.3 km. 시험 비행에서 19 km쯤에 경고가 뜬 것과 맞습니다. 14 m/s 근처에서는 23.6 km까지 갑니다. 배터리를 바꾸지 않아도 섬에 갈 수 있다는 뜻입니다.

꼭대기를 더 정확히 찾는 가장 단순한 방법은 아주 촘촘하게 다 계산해 보는 격자 탐색(grid search) 입니다. 0.01 m/s 간격이면 점이 2001개입니다.

>>> vs = np.linspace(5, 25, 2001)
>>> R = range_km(vs)
>>> round(vs[np.argmax(R)], 2), round(R.max(), 2)
(np.float64(14.12), np.float64(23.57))
>>> round(vs[np.argmax(flight_min(vs))], 2)
np.float64(10.17)
>>> round(flight_min(14.12), 1), round(range_km(10.17), 1)
(27.8, 19.8)

np.argmax는 가장 큰 값의 위치(번호)를 돌려줍니다. 여기서 반전이 하나 있습니다. 가장 오래 떠 있는 속도(10.17 m/s, 32.4분)와 가장 멀리 가는 속도(14.12 m/s, 23.6 km)가 다릅니다. 오래 떠 있기만 하면 되는 감시 비행이라면 10 m/s, 멀리 가야 하는 배송이라면 14 m/s입니다. 무엇을 최대로 할지 먼저 정해야 합니다.

격자 탐색은 확실하지만 비쌉니다. 0.01 m/s 정확도에 2001번 계산했으니, $10^{-6}$까지 가려면 2000만 번입니다. 이번 함수는 한 번 계산이 순식간이지만, 실무에서는 한 번 계산이 비행 시뮬레이션 한 번(몇 초~몇 분)인 경우가 흔합니다. 그래서 함수를 몇 번 계산했는가가 최적화 방법의 성적표가 됩니다.

6-2. 꼭대기는 기울기가 0인 곳: 근 찾기로 바꾸기 (7분)

언덕 꼭대기에 서면 앞으로도 뒤로도 평평합니다. 매끄러운 함수의 꼭대기에서는 기울기가 0입니다.

그러면 최적화는 4·5교시의 근 찾기 문제가 됩니다. "$R'(v) = 0$인 $v$를 찾아라." 도함수 식을 손으로 구하기 번거로우니 3교시의 중앙 차분으로 근사하고, 5교시의 뉴턴 방법을 씁니다. 뉴턴 방법은 $g(v) = 0$을 풀 때 $g'(v)$가 필요하니, 여기서는 $R''(v)$(2계 도함수)까지 차분으로 근사합니다.

>>> def d1(f, x, h=1e-3):
...     return (f(x + h) - f(x - h)) / (2 * h)
>>> def d2(f, x, h=1e-3):
...     return (f(x + h) - 2 * f(x) + f(x - h)) / h**2
>>> v = 10.0
>>> for k in range(5):
...     v = v - d1(range_km, v) / d2(range_km, v)
...     print(k + 1, round(v, 6))
1 14.589979
2 14.08973
3 14.115919
4 14.115988
5 14.115988
>>> v_best = v
>>> round(d1(range_km, v_best), 9), round(d2(range_km, v_best), 4)
(0.0, -0.4114)

네 번 만에 14.115988에 도착했습니다. 맞는 자릿수가 2개, 4개, 7개로 늘어나는 5교시의 2차 수렴 그대로입니다. 마지막 줄의 $R'' < 0$은 이곳이 위로 볼록한 꼭대기라는 확인입니다.

이 방법에는 함정이 둘 있습니다. 첫째, 한 걸음에 함수를 세 번 계산하고 도함수를 두 번이나 근사합니다. 둘째, 기울기 0은 꼭대기만의 특징이 아닙니다. 골짜기 바닥도 기울기가 0입니다. 봉우리가 둘인 함수에서 두 봉우리 사이에서 출발해 봅시다.

>>> def bumps(x):
...     return 2.5 * np.exp(-(x - 2.2)**2 / 4) + 3 * np.exp(-(x - 8)**2 / 1.2)
>>> x = 5.0
>>> for k in range(6):
...     x = x - d1(bumps, x) / d2(bumps, x)
>>> round(x, 4), round(d2(bumps, x), 4)
(np.float64(5.7788), np.float64(0.8666))

뉴턴 방법은 성실하게 기울기 0인 곳을 찾았지만, 2계 도함수가 양수입니다. 골짜기 바닥(최솟값) 에 도착한 것입니다. 뉴턴 방법은 꼭대기와 골짜기를 구별하지 않습니다. 그래서 도함수 없이, 함숫값의 크기 비교만으로 꼭대기를 쫓는 방법이 필요합니다.

6-3. 단봉 함수와 두 점 비교: 이분법이 바로 안 되는 이유 (6분)

4교시의 이분법은 가운데 점 하나로 구간을 반으로 줄였습니다. 근이 있는 쪽은 부호로 알 수 있었습니다. 최댓값에서는 이 방법이 통하지 않습니다. 가운데 점의 값 $f(m)$ 하나만으로는 꼭대기가 왼쪽에 있는지 오른쪽에 있는지 알 수 없기 때문입니다. 값이 크다는 것만 알지, 오르막인지 내리막인지는 모릅니다.

그래서 구간 안에 점을 두 개 찍습니다. 단, 구간 $[a, b]$ 안에 봉우리가 하나뿐이어야 합니다. 꼭대기까지는 계속 오르고, 꼭대기를 지나면 계속 내려가는 함수를 단봉 함수(unimodal function) 라고 합니다. 단봉 함수에서 $a < x_1 < x_2 < b$인 두 점을 비교하면 이렇게 판단할 수 있습니다.

가장 먼저 떠오르는 방법은 구간을 3등분하는 두 점을 쓰는 것입니다. 이것을 3분 탐색(ternary search) 이라고 합니다.

>>> def ternary_max(f, a, b, n):
...     for k in range(n):
...         x1, x2 = a + (b - a) / 3, b - (b - a) / 3
...         if f(x1) > f(x2):
...             b = x2
...         else:
...             a = x1
...     return a, b
>>> a, b = ternary_max(range_km, 5, 25, 10)
>>> round(a, 4), round(b, 4), round(b - a, 4)
(13.9356, 14.2825, 0.3468)
>>> round(20 * (2 / 3)**10, 4)
0.3468

구간이 단계마다 2/3로 줄어 10단계 만에 폭 0.35가 되었습니다. 그런데 단계마다 새 점 두 개를 계산하니 함수를 20번 계산했습니다. 계산 한 번당 구간은 $\sqrt{2/3} \approx 0.816$배밖에 줄지 않습니다. 남은 구간 안에는 지난 단계에서 계산한 점 하나가 남아 있는데, 그 값을 버리고 새로 두 점을 찍기 때문입니다. 그 점을 다시 쓸 수 있게 두 점의 위치를 고르면 어떨까요? 다음 절의 황금분할 탐색이 바로 그 생각입니다.

6-4. 황금분할 탐색: 한 점을 다시 쓰는 비율 (12분)

구간 $[a, b]$의 길이를 $L$이라 하고, 두 안쪽 점을 양 끝에서 같은 비율 $r$만큼 떨어진 곳에 둡니다.

$f(x_1) > f(x_2)$라서 $(x_2, b]$를 버렸다고 합시다. 새 구간은 $[a, x_2]$, 길이는 $rL$입니다. 이 새 구간의 오른쪽 안쪽 점은 $a + r \cdot rL$에 있어야 합니다. 그런데 이 자리가 마침 옛 $x_1 = a + (1 - r)L$과 같다면, 옛 $x_1$의 값을 그대로 다시 쓸 수 있습니다. 조건을 적으면 이렇습니다.

이 비율이 고대 그리스부터 알려진 황금비(golden ratio) $\varphi \approx 1.618$의 역수 $1/\varphi = 0.618\ldots$입니다. 그래서 이 방법을 황금분할 탐색(golden-section search) 이라고 부릅니다.

>>> r = (np.sqrt(5) - 1) / 2
>>> round(r, 6), round(r**2, 6), round(1 - r, 6)
(np.float64(0.618034), np.float64(0.381966), np.float64(0.381966))

$r^2$과 $1 - r$이 정말 같습니다. 정리하면 황금분할 탐색은 단계마다 함수를 한 번만 새로 계산하고, 구간을 0.618배로 줄입니다. 3분 탐색의 계산 한 번당 0.816배보다 훨씬 빠릅니다.

손으로 첫 단계를 해 봅시다. $[5, 25]$에서 $x_1 = 25 - 0.618 \times 20 = 12.64$, $x_2 = 5 + 0.618 \times 20 = 17.36$입니다. $R(12.64) \approx 23.0$, $R(17.36) \approx 21.9$이므로 $f(x_1) > f(x_2)$, 새 구간은 $[5, 17.36]$입니다. 옛 $x_1 = 12.64$는 새 구간의 $x_2$가 되고, 새로 계산할 점은 $x_1 = 17.36 - 0.618 \times 12.36 = 9.72$ 하나뿐입니다.

이제 파이썬으로 짜서 표를 찍어 봅니다. calls는 함수 계산 횟수입니다. 처음 두 점에 2번, 그 뒤로 단계마다 1번입니다.

>>> def golden_max(f, a, b, tol=1e-3, show=False):
...     x1, x2 = b - r * (b - a), a + r * (b - a)
...     f1, f2 = f(x1), f(x2)
...     calls, k = 2, 0
...     while b - a > tol:
...         if f1 > f2:
...             b, x2, f2 = x2, x1, f1
...             x1 = b - r * (b - a)
...             f1 = f(x1)
...         else:
...             a, x1, f1 = x1, x2, f2
...             x2 = a + r * (b - a)
...             f2 = f(x2)
...         calls, k = calls + 1, k + 1
...         if show:
...             print(f"{k:2d} {a:8.4f} {b:8.4f} {x1:8.4f} {x2:8.4f} {b - a:7.4f}")
...     return (a + b) / 2, calls
>>> x, calls = golden_max(range_km, 5, 25, tol=0.01, show=True)
 1   5.0000  17.3607   9.7214  12.6393 12.3607
 2   9.7214  17.3607  12.6393  14.4427  7.6393
 3  12.6393  17.3607  14.4427  15.5573  4.7214
 4  12.6393  15.5573  13.7539  14.4427  2.9180
 5  13.7539  15.5573  14.4427  14.8684  1.8034
 6  13.7539  14.8684  14.1796  14.4427  1.1146
 7  13.7539  14.4427  14.0170  14.1796  0.6888
 8  14.0170  14.4427  14.1796  14.2801  0.4257
 9  14.0170  14.2801  14.1175  14.1796  0.2631
10  14.0170  14.1796  14.0791  14.1175  0.1626
11  14.0791  14.1796  14.1175  14.1412  0.1005
12  14.0791  14.1412  14.1028  14.1175  0.0621
13  14.1028  14.1412  14.1175  14.1266  0.0384
14  14.1028  14.1266  14.1119  14.1175  0.0237
15  14.1119  14.1266  14.1175  14.1210  0.0147
16  14.1119  14.1210  14.1154  14.1175  0.0091
>>> round(x, 4), calls
(np.float64(14.1164), 18)

표에서 점이 다시 쓰이는 모습을 확인해 보세요. 1단계의 $x_2 = 12.6393$이 2단계에서는 $x_1$이 되고, 2단계의 $x_2 = 14.4427$은 3단계의 $x_1$이 됩니다. 마지막 열의 폭은 매번 정확히 0.618배입니다. 함수 18번 계산으로 격자 탐색(2001번)과 같은 정확도에 도착했습니다.

구간이 매번 $r$배가 되니, 처음 길이 $L$을 허용 오차 $\text{tol}$ 아래로 줄이는 데 필요한 단계 수 $n$은 미리 계산할 수 있습니다.

>>> L, tol = 20, 0.01
>>> int(np.ceil(np.log(tol / L) / np.log(r)))
16
>>> int(np.ceil(np.log(1e-6 / L) / np.log(r)))
35

위의 표와 같은 16단계입니다. 정확도를 $10^{-6}$까지 올려도 35단계, 함수 계산 37번이면 됩니다. 4교시의 이분법처럼 "몇 번 돌리면 되는지"를 돌리기 전에 아는 방법입니다. 다만 이분법은 단계마다 1/2로 줄었으니, 황금분할 탐색은 이분법보다 조금 느립니다. 부호 같은 방향 정보 없이 크기 비교만 쓰는 대가입니다.

아래 실습대의 남색 곡선이 항속 거리 R(v)이고, 회색 점선이 참 꼭대기입니다. 옅은 남색 띠가 지금 구간 [a, b]이고, 구간 밖은 회색으로 덮입니다. 호박색 점이 x₁, 청록 점이 x₂입니다. '다음'이나 '재생'을 누르면 한 단계씩 나아가고, 방금 버린 구간이 붉게 깜박였다가 사라집니다. 그래프 바닥의 두 손잡이 a₀, b₀를 끌면(또는 화살표 키) 처음 구간이 바뀝니다. 아래 표에는 단계마다 a, b, x₁, x₂, 두 점의 값, 구간 폭이 쌓입니다. 구간이 너무 좁아져 안 보이면 '지금 구간 확대해 보기'를 켜세요.

6-5. 포물선 보간: 세 점으로 꼭대기 짐작하기 (6분)

황금분할 탐색은 함숫값의 크기 비교만 씁니다. $R(x_1) = 23.0$이든 $23.5$든 "더 크다"는 사실만 남기고 숫자는 버립니다. 숫자까지 쓰면 더 빨리 갈 수 있지 않을까요?

매끄러운 함수의 꼭대기 근처는 포물선(2차 함수)과 많이 닮았습니다. 그렇다면 세 점 $a < m < b$를 지나는 포물선을 그리고, 그 꼭짓점으로 바로 뛰어가 봅시다. 이것이 포물선 보간(parabolic interpolation) 입니다. 세 점을 지나는 포물선의 꼭짓점은 다음 식으로 구합니다.

분모가 양수이면 포물선이 위로 볼록(∩)해서 꼭짓점이 꼭대기입니다. 분모가 0 이하이면 아래로 볼록(∪)이거나 일직선이라 꼭짓점이 꼭대기가 아닙니다. 그때는 멈춥니다. 꼭짓점 $x^*$의 값을 계산한 뒤에는 꼭대기를 감싸는 세 점을 다시 고릅니다($x^*$와 $m$ 중 큰 쪽이 새 가운데 점).

>>> def parabola_max(f, a, b, tol=1e-6, max_iter=30):
...     m = (a + b) / 2
...     fa, fm, fb = f(a), f(m), f(b)
...     calls = 3
...     for k in range(max_iter):
...         num = (m - a)**2 * (fm - fb) - (m - b)**2 * (fm - fa)
...         den = (m - a) * (fm - fb) - (m - b) * (fm - fa)
...         if den <= 0:
...             return None, calls
...         x = m - 0.5 * num / den
...         if not a < x < b:
...             return None, calls
...         fx = f(x)
...         calls += 1
...         if abs(x - m) < tol:
...             return x, calls
...         if x > m:
...             if fx > fm:
...                 a, fa, m, fm = m, fm, x, fx
...             else:
...                 b, fb = x, fx
...         else:
...             if fx > fm:
...                 b, fb, m, fm = m, fm, x, fx
...             else:
...                 a, fa = x, fx
...     return m, calls
>>> xp, cp = parabola_max(range_km, 5, 25)
>>> round(xp, 6), cp, f"{abs(xp - v_best):.1e}"
(14.115988, 10, '1.0e-07')
>>> xg, cg = golden_max(range_km, 5, 25, tol=1e-6)
>>> round(xg, 6), cg, f"{abs(xg - v_best):.1e}"
(np.float64(14.115988), 37, '9.3e-08')

같은 정확도에 포물선 보간은 함수 계산 10번, 황금분할 탐색은 37번입니다. 꼭대기 근처에서 함수가 포물선을 닮을수록 포물선 보간은 무섭게 빨라집니다. 그런데 봉우리가 두 개인 함수에서는 어떨까요?

>>> parabola_max(bumps, 0, 10)
(None, 3)
>>> golden_max(bumps, 0, 10)
(np.float64(2.199876131282677), 22)
>>> golden_max(bumps, 4, 10)
(np.float64(7.999546896146215), 21)

포물선 보간은 첫 세 점(0, 5, 10)이 아래로 볼록한 포물선을 만들어 첫 걸음에 멈췄습니다. 꼭짓점이 구간 밖으로 튀어 나가는 경우도 있습니다. 황금분할 탐색은 늘 무언가를 찾지만, $[0, 10]$에서는 낮은 봉우리(x = 2.2, 높이 2.5)를 찾았습니다. 구간 안이 단봉이라는 약속이 깨졌기 때문입니다. 구간을 $[4, 10]$으로 잡아야 높은 봉우리(x = 8, 높이 3)를 찾습니다. 이렇게 주변보다는 높지만 전체에서 가장 높지는 않은 곳을 국소 최댓값(local maximum), 전체에서 가장 높은 곳을 전역 최댓값(global maximum) 이라고 합니다.

위 실습대에서 방법을 '포물선 보간'으로 바꿔 보세요. 호박색 세 점을 지나는 청록 점선 포물선이 그려지고, 그 꼭짓점 x*로 뛰어갑니다. '봉우리 두 개'에서는 포물선이 아래로 볼록해지는 모습도 볼 수 있습니다. 오른쪽 아래 줄에 두 방법이 참 꼭대기에 0.001 이내로 다가가기까지의 계산 횟수가 나옵니다.

6-6. 변수가 둘이면: 경사하강법 맛보기 (5분)

순항 고도도 정할 수 있다면 어떨까요? 낮게 날면 지형 때문에 난기류가 많고, 높이 날면 맞바람이 셉니다. 관제실이 지난달 비행 기록으로 1 km당 쓰는 에너지 $e(v, h)$(Wh/km)를 근사식으로 맞춰 두었다고 합시다. 고도 $h$는 100 m 단위입니다. 이런 식을 데이터로 맞추는 방법은 9교시에서 배웁니다.

1 km당 에너지가 가장 적은 곳이 가장 멀리 가는 곳입니다(거리 = $E / e$). 이번에는 최소화입니다. 변수가 둘이면 "구간"을 좁힐 수 없습니다. 대신 지금 자리에서 가장 가파르게 내려가는 방향으로 조금씩 걸어갑니다. 가장 가파르게 올라가는 방향이 그레이디언트(gradient) $\nabla e = (\partial e/\partial v,\ \partial e/\partial h)$이므로, 그 반대 방향으로 갑니다.

이것이 경사하강법(gradient descent) 입니다. $\eta$(에타)는 한 걸음의 크기를 정하는 학습률(learning rate) 입니다. 기울기는 변수마다 중앙 차분으로 구합니다(3교시).

>>> def energy(p):
...     v, h = p
...     return 17 + 0.15 * (v - 14)**2 + 0.8 * (h - 0.9)**2 + 0.1 * (v - 14) * (h - 0.9)
>>> def grad(f, p, h=1e-5):
...     g = np.zeros(len(p))
...     for i in range(len(p)):
...         e = np.zeros(len(p))
...         e[i] = h
...         g[i] = (f(p + e) - f(p - e)) / (2 * h)
...     return g
>>> p = np.array([20.0, 0.3])
>>> np.round(grad(energy, p), 4)
array([ 1.74, -0.36])
>>> for k in range(3):
...     p = p - 1.0 * grad(energy, p)
...     print(k + 1, np.round(p, 3), round(energy(p), 3))
1 [18.26  0.66] 19.666
2 [17.006  0.618] 18.334
3 [16.132  0.769] 17.668

지금 운항 조건(20 m/s, 고도 30 m)에서 기울기는 (1.74, −0.36)입니다. 속도는 줄이고 고도는 올리라는 뜻입니다. 학습률을 바꿔 30걸음씩 가 봅시다.

>>> def descend(eta, steps=30):
...     p = np.array([20.0, 0.3])
...     for _ in range(steps):
...         p = p - eta * grad(energy, p)
...     return np.round(p, 3), round(energy(p), 3)
>>> descend(0.3)
(array([14.383,  0.871]), np.float64(17.022))
>>> descend(1.0)
(array([14. ,  0.9]), np.float64(17.0))
>>> descend(1.5)
(array([ -317.779, -4337.604]), np.float64(15218565.383))

$\eta = 0.3$은 안전하지만 30걸음으로는 아직 덜 왔습니다. $\eta = 1.0$은 (14 m/s, 90 m)에 정확히 도착했습니다. $\eta = 1.5$는 골짜기 건너편을 넘나들며 점점 크게 튀어 나가, 30걸음 만에 숫자가 수천 단위가 되었습니다(고도 −4337은 땅속 43만 m입니다). 1교시에서 걸음이 너무 큰 오일러 방법이 출렁이다 튀어 나간 것과 같은 현상입니다. 학습률이 작으면 느리고, 크면 발산합니다. 기계 학습에서 학습률을 조심스럽게 고르는 이유가 바로 이것입니다.

책의 기호 읽는 법

실무에서는

SciPy의 scipy.optimize가 이번 교시의 도구를 모두 갖추고 있습니다. 모두 최소화 함수이니 최댓값은 부호를 뒤집어 넣습니다.

from scipy.optimize import minimize_scalar, minimize

# 1차원: 구간을 주면 브렌트 방법 (황금분할 + 포물선 보간)
res = minimize_scalar(lambda v: -range_km(v), bounds=(5, 25), method="bounded")
res.x, -res.fun, res.nfev      # 꼭대기 위치, 최댓값, 함수 계산 횟수

# 황금분할만 쓰고 싶을 때 (bracket은 출발 구간일 뿐, 그 밖으로 나갈 수도 있다)
minimize_scalar(lambda v: -range_km(v), bracket=(5, 25), method="golden")

# 변수가 여럿: 기본은 BFGS (기울기로 방향과 걸음을 함께 고르는 경사하강법의 발전형)
res = minimize(energy, x0=[20, 0.3])
res.x

minimize_scalar의 기본 방법인 브렌트 방법(Brent's method) 은 이번 교시의 두 방법을 섞은 것입니다. 포물선 보간으로 빠르게 가다가, 포물선이 길을 잃을 것 같으면 황금분할 탐색 한 걸음으로 물러납니다. 5교시에서 근 찾기용 브렌트 방법(brentq)이 이분법과 할선법을 섞었던 것과 같은 생각입니다. 봉우리가 여럿일 수 있으면 먼저 np.argmax로 격자를 훑어 큰 그림을 보고, 그 근처에서 이 함수들을 부르세요.

자주 하는 실수 (3분)

최소화 함수에 부호를 안 뒤집고 넣기. minimize_scalar(range_km, ...)는 가장 짧게 가는 속도를 찾아 구간 끝(5 m/s)을 돌려줍니다. 최댓값은 lambda v: -range_km(v)로 넣으세요.

찾은 값을 전역 최댓값이라고 믿기. 황금분할 탐색과 포물선 보간은 구간 안이 단봉이라고 믿고 움직입니다. 봉우리가 여럿이면 국소 최댓값에 멈춥니다. 먼저 격자로 훑어 그래프를 보세요.

꼭대기의 x를 지나치게 정밀하게 원하기. 꼭대기 근처는 평평해서, x가 꽤 달라도 함숫값은 거의 같습니다.

>>> round(range_km(v_best + 0.1) - range_km(v_best), 4)
-0.002
>>> range_km(v_best + 1e-8) - range_km(v_best)
3.552713678800501e-15

순항 속도를 0.1 m/s 틀려도 거리 차이는 2 m입니다. 그래서 관제실은 14 m/s로 반올림해 써도 됩니다. 반대로 $10^{-8}$만큼 옮기면 차이가 반올림 오차(2교시)만 한 크기입니다. 꼭대기에서 벗어났는데 오히려 값이 커졌다고 나옵니다. 컴퓨터가 어느 쪽이 높은지 구별하지 못하는 것입니다. 그래서 꼭대기의 위치는 대략 $\sqrt{\varepsilon} \approx 10^{-8}$의 상대 정확도 아래로는 찾을 수 없습니다. 허용 오차를 그보다 작게 잡으면 계산만 늘어납니다.

연습 문제 (3분)

1. 짐을 더 실어 유도 전력 계수 A가 4500에서 6000이 되었습니다. golden_max로 가장 멀리 가는 속도와 그 거리를 구해 보세요. 속도는 빨라질까요, 느려질까요?
2. 감시 비행처럼 오래 떠 있는 것이 목표라면, golden_max로 비행 시간 flight_min의 꼭대기를 허용 오차 1e-4로 찾아 보세요. 함수 계산 횟수는 단계 수 공식과 맞나요?
3. 1 km당 에너지 power(v) / (3.6 * v) (Wh/km)를 최소로 하는 속도를 구하세요. golden_max는 최댓값만 찾으니 부호를 뒤집어 넣습니다. 답이 항속 거리의 꼭대기와 같은지 확인해 보세요.

풀이

>>> def range_heavy(v):
...     return 3.6 * v * E / (6000 / v + B * v**3 + C)
>>> x, calls = golden_max(range_heavy, 5, 25, tol=1e-4)
>>> round(x, 3), round(range_heavy(x), 2), calls
(np.float64(15.021), np.float64(21.12), 28)
>>> x, calls = golden_max(flight_min, 5, 25, tol=1e-4)
>>> round(x, 3), round(flight_min(x), 1), calls
(np.float64(10.174), np.float64(32.4), 28)
>>> int(np.ceil(np.log(1e-4 / 20) / np.log(r))) + 2
28
>>> x, calls = golden_max(lambda v: -power(v) / (3.6 * v), 5, 25, tol=1e-4)
>>> round(x, 3), round(power(x) / (3.6 * x), 3), round(E / (power(x) / (3.6 * x)), 2)
(np.float64(14.116), np.float64(16.974), np.float64(23.57))

1번은 가장 좋은 속도가 14.12에서 15.02 m/s로 빨라지고, 거리는 23.57에서 21.12 km로 줄어듭니다. 무거워지면 몸을 띄우는 전력(A/v)의 몫이 커지는데, 이 항은 빨리 날수록 작아지기 때문입니다. 섬 보건소(21 km)까지는 아슬아슬합니다.

2번은 단계 수 공식으로 26단계에 처음 두 점 2번을 더해 28번입니다.

3번은 1 km당 16.974 Wh를 쓰는 14.116 m/s가 답입니다. 1 km당 에너지 $e(v) = P(v)/(3.6v)$를 최소로 하는 속도가 $R(v) = E/e(v)$를 최대로 하는 속도와 같습니다. $E$로 나눈 값의 최대는 나누는 값의 최소이기 때문입니다. 배터리 용량 $E$가 답에 들어가지 않는다는 점도 눈여겨보세요. 배터리를 바꿔도 가장 멀리 가는 속도는 그대로입니다.

다음 날 아침, 순항 속도를 14 m/s로 바꾼 3호기가 섬 보건소에 무사히 도착했습니다. 그런데 민재 팀장의 얼굴이 밝지 않습니다. "속도를 늦추니까 하루에 나를 수 있는 물량이 바뀌었어요. 항구 허브, 읍내 허브, 산마을 허브 세 곳의 들고 나는 물량이 서로 맞게 배정량을 다시 짜야 해요." 다음 교시에는 미지수가 여러 개인 방정식을 한꺼번에 푸는 가우스 소거를 배웁니다.