제1부 수학적 기초

1.13 최적화 II: 볼록 최적화와 QP — 활성 집합·ADMM·내점법·qpOASES·OSQP·ProxQP

1.13.1 이 장에서 배우는 것 / 선수 지식 / 특히 중요한 독자

이 장은 선형 제약 아래에서 볼록 이차 함수를 최소화하는 이차 계획법(quadratic programming, QP)과 그 해법기를 다룬다. 미분 역기구학의 관절 속도 한계, 접촉력 배분의 마찰 제약, 선형 MPC의 입력·상태 제약은 모두 QP로 적히며, 제어기는 이 QP를 제어 주기마다 푼다. 앞 장(1.12)이 제약 없는 비선형 최소제곱을 다루었다면, 이 장은 제약을 포함하는 볼록 문제를 다룬다. 다룰 질문은 다섯 가지다. QP의 최적성은 어떤 조건으로 판정하는가, 활성 집합법·내점법·ADMM은 각각 무엇을 반복하는가, qpOASES·OSQP·ProxQP는 설계의 어느 부분이 다른가, 열 시작(warm start, 이전 QP의 해로 시작)과 $\rho$·스케일링은 반복 수를 얼마나 바꾸는가, 해법기가 반환하는 '해결됨' 상태를 얼마나 신뢰할 수 있는가.

배우는 것:

• 볼록 QP와 KKT — $\ell\le Ax\le u$ 형태, 필요충분조건, 승수의 해석, 쌍대성 간격.
• 활성 집합법 — 작업 집합, 등식 QP의 KKT 선형계, Goldfarb–Idnani 쌍대법, qpOASES의 온라인 활성 집합.
• 내점법 — 중심 경로, 원-쌍대 뉴턴 단계, Mehrotra 예측-수정.
• ADMM(OSQP) — 변수 분할, $\rho\cdot\sigma\cdot\alpha$, 정지 기준, 스케일링·폴리싱·열 시작.
• 근접 증강 라그랑주(ProxQP) — 외부 근접 곱셈자 반복과 내부 준매끄러운 뉴턴 반복.
• 측정 — 직접 구현한 세 해법의 KKT 잔차·$\rho$ 민감도·열 시작과 OSQP·ProxQP·qpOASES 등 실제 라이브러리와의 대조.

선수 지식: 1.1(행렬 분해), 1.2(숄레스키·조건수·감쇠 최소제곱), 1.12. 표기는 BK-05와 A.2을 따른다.

특히 중요한 독자: QP 해법기를 선택하고 설정하는 제어·소프트웨어 개발자, 전신 제어·MPC 연구자, 안전 기능에 QP를 사용하는 시스템의 인증 담당 엔지니어(상태 코드의 해석).

출처 방침. 서지는 Crossref 응답으로 확인했다(목록은 이 장 끝). 직접 읽은 자료: ① OSQP C 소스(osqp/osqp master — 기본값 상수·ADMM 갱신·정지 판정·$\rho$ 갱신·폴리싱), ② qpOASES 3.2 사용자 매뉴얼(2017-04)과 Options.cpp, ③ ProxQP 논문(RSS 2022, HAL 판)과 ProxSuite의 settings.hpp·solver.hpp·doc/*.md, ④ CasADi의 qpoases_interface.cpp, ⑤ CVXPY의 osqp_qpif.py. 직접 실행한 것: osqp 1.1.3, proxsuite 0.8.2, casadi 3.8.2(번들 qpOASES), quadprog 0.1.14, clarabel 0.11.1, cvxpy 1.9.3, scipy 1.18.1(numpy 2.5.4, 격리된 venv). 읽지 못한 것: 인용한 교과서·논문의 본문(Goldfarb–Idnani, Mehrotra, OSQP 논문, ADMM 수렴 이론 — 정리·장 번호는 (확인 필요)), 상용 해법기. 이 장의 수치는 모두 순수 파이썬 구현으로 계산했으며 1.1의 LU 풀이와 1.2의 숄레스키를 불러와 사용했다.

1.13.2 개념 — 제약이 있는 최소제곱

1.2의 최소제곱은 $\min\tfrac12\lVert Jx-b\rVert^2$이며 제약이 없다. 관절 속도 한계 $\lvert\dot q_i\rvert\le\dot q_{\max}$를 더하면 목적은 이차 그대로이고 제약만 선형이므로 QP가 된다. 제약이 있으면 해는 일반적으로 닫힌 식으로 주어지지 않는다. 최적점에서 일부 제약은 등호로 성립하고(활성), 나머지는 해에 영향을 주지 않는다. 해법은 활성 집합을 찾는 방식에 따라 세 계열로 나뉜다.

1.13.3 원리

1.13.3.1 볼록 QP의 형태

결론: 변수 경계·등식·일방 부등식은 모두 $\ell\le Ax\le u$ 한 형태로 적을 수 있고, $P\succeq0$이면 모든 국소 최적점이 전역 최적점이다.

$$\min_{x\in\mathbb R^n}\ \tfrac12x^\top Px+q^\top x\quad\text{s.t.}\quad\ell\le Ax\le u,\qquad P\succeq0$$

등식은 $\ell_i=u_i$, 일방 제약은 $\ell_i=-\infty$ 또는 $u_i=+\infty$, 변수 경계는 $A$의 단위행렬 행으로 적는다. 이 형태는 OSQP의 입력이다. qpOASES는 $\tfrac12x^\top Hx+x^\top g$와 $lbA\le Ax\le ubA$·$lb\le x\le ub$(매뉴얼 1.1절), ProxQP는 $Ax=b$와 $l\le Cx\le u$(ProxSuite README)를 받는다. 이 장은 OSQP 형태로 통일한다. $P\succ0$이면 해가 유일하다. $P=0$이면 선형 계획법이고, 제약이 없고 $P=J^\top J+\lambda I$이면 감쇠 최소제곱(1.2.3.5)이다. QP를 포함하는 더 넓은 볼록 문제 계층(2차 원뿔·반정부호 계획법)은 1.17가, $P$가 불정부호인 비볼록 QP는 1.14가 다룬다.

1.13.3.2 KKT 조건과 쌍대성

결론: 선형 제약 볼록 QP에서 KKT 조건은 최적성의 필요충분조건이고 승수 $y$는 제약 경계에 대한 최적값의 민감도이다. 따라서 해법기의 상태 코드와 별개로 KKT 잔차를 직접 계산해 해를 판정할 수 있다.

$$Px+q+A^\top y=0,\qquad Ax=z,\qquad z\in[\ell,u],\qquad y_i\ \begin{cases}\ge0&z_i=u_i\\ \le0&z_i=\ell_i\\ =0&\ell_i<z_i<u_i\end{cases}$$

ProxQP 논문은 선형 제약 볼록 QP에서 강쌍대성이 성립하고 KKT 조건이 필요충분임을 서술한다. 이 장의 부호는 OSQP를 따른다(상한이 활성이면 $y_i>0$). 판정에 사용하는 잔차는 원 잔차 $r_{\rm prim}=\lVert Ax-z\rVert_\infty$, 쌍대 잔차 $r_{\rm dual}=\lVert Px+q+A^\top y\rVert_\infty$, 쌍대성 간격이다.

$$g=x^\top Px+q^\top x+S(y),\qquad S(y)=u^\top\max(y,0)+\ell^\top\min(y,0)$$

$g$는 원 목적과 쌍대 목적의 차이며 상보성을 포함한다. 비퇴화 문제에서 승수는 $\partial f^\ast/\partial u_i=-y_i$로 읽는다. 접촉력 배분에서 마찰 제약의 승수가 0이 아니면 그 접촉이 마찰 한계에 있다. $P\succ0$이고 제약이 $Cx\le d$이면 쌍대 문제는 $\max_{\lambda\ge0}\,-\tfrac12(q+C^\top\lambda)^\top P^{-1}(q+C^\top\lambda)-d^\top\lambda$이며, 쌍대 활성 집합법은 이 문제 위에서 진행한다 (확인 필요: Goldfarb–Idnani 1983의 정리).

1.13.3.3 활성 집합법

결론: 활성 집합법은 작업 집합 $\mathcal W$를 등식으로 둔 QP의 KKT 선형계 하나를 풀고, 가능성 또는 승수 부호를 위반한 제약 하나를 작업 집합에 더하거나 뺀다. 퇴화하지 않는 문제에서는 유한 반복에 기계 정밀도의 정확한 해를 준다.

제약을 $a_i^\top x\ge b_i$로 적을 때 작업 집합의 등식 QP는 다음 선형계이다.

$$\begin{bmatrix}P&A_{\mathcal W}^\top\\A_{\mathcal W}&0\end{bmatrix}\begin{bmatrix}x\\-\mu\end{bmatrix}=\begin{bmatrix}-q\\b_{\mathcal W}\end{bmatrix},\qquad \mu\ge0$$

• 주 활성 집합법 — 가능한 점에서 시작해 차단 제약을 더하고 $\mu_j<0$인 제약을 뺀다. 가능한 시작점이 필요하다. ProxQP 논문은 입문 문헌으로 Nocedal & Wright의 16.5절을 가리킨다 (확인 필요: 본문).
• 쌍대 활성 집합법(Goldfarb–Idnani, 1983) — $P\succ0$을 가정하고 제약 없는 최소점에서 출발해 가장 크게 위반된 제약을 하나씩 더한다. 쌍대 가능성 $\mu\ge0$을 유지하며 원 가능성을 얻는다. 가능한 시작점이 필요 없고 불가능한 문제는 단계 길이가 무한대가 되어 검출된다. quadprog가 이 알고리즘의 구현이다(PyPI 메타데이터).
• 매개변수 활성 집합법(qpOASES) — 기울기·경계가 매개변수에 아핀으로 의존할 때 직전 QP의 최적해에서 새 QP로 호모토피를 따라 이동한다(매뉴얼 1.1·5.4절). 활성 집합이 바뀌는 점마다 작업 집합이 재계산되며 nWSR가 이 횟수이다. 연속한 QP에서 최적 활성 집합이 크게 변하지 않는다는 기대에 기초한다(매뉴얼 1.1절).

단점은 문헌과 매뉴얼에서 읽은 범위에서 세 가지다. 활성 집합의 순환(cycling)과 제약 규약 위배 시 둔화(ProxQP 논문 2절), 밀집 선형대수 중심의 구현(qpOASES 매뉴얼 4.7절은 밀집 행렬을 가진 중소 규모 QP를 위해 개발되었다고 밝힌다), 반복 수의 최악 경계가 문제 크기에 대해 지수적일 수 있다는 알려진 결과 (확인 필요: 출처)이다.

1.13.3.4 원-쌍대 내점법

결론: 내점법은 부등식을 슬랙 변수와 상보성 섭동으로 처리하여 가능 영역의 내부를 지나며, 반복마다 KKT 선형계 하나를 풀고, 반복 수가 문제 크기에 거의 의존하지 않는다. ProxQP 논문 2절은 호모토피 구조 때문에 열 시작에 사용하기 어렵다는 점을 내점법의 주된 단점으로 든다.

$Gx+s=h$, $Ex=f$, $s,\lambda>0$인 형태에서 섭동된 KKT 조건은 다음과 같다. $\mu=s^\top\lambda/m$은 상보 간격의 평균이고 $\sigma\in[0,1]$은 중심화 정도이다.

$$Px+q+G^\top\lambda+E^\top\nu=0,\qquad Ex=f,\qquad Gx+s=h,\qquad s\circ\lambda=\sigma\mu\mathbf 1$$

슬랙 $s$와 승수 $\lambda$의 증분을 소거하면 뉴턴 단계는 $D=\operatorname{diag}(\lambda_i/s_i)$에 대한 행렬 $\begin{bmatrix}P+G^\top DG&E^\top\\E&0\end{bmatrix}$의 풀이로 귀착된다. $D$가 반복마다 변하므로 분해를 매번 다시 한다. $\sigma=1$이면 뉴턴 단계는 현재 $\mu$의 중심 경로 점을 향한다. 중심 경로는 로그 장벽 문제 $\min\ \tfrac12x^\top Px+q^\top x-\mu\sum_i\log(h_i-g_i^\top x)$의 해 $x(\mu)$이며 $\mu\to0$에서 QP의 해에 수렴한다. Mehrotra의 예측-수정은 $\sigma=0$인 예측 단계의 결과로 $\sigma=(\mu_{\rm aff}/\mu)^3$을 정하고 2차 항을 포함한 수정 단계를 한 번 더 푼다 (확인 필요: 원 논문의 식).

1.13.3.5 ADMM — OSQP의 분할·$\rho$·정지 기준

결론: OSQP는 QP를 $Ax=z$, $z\in[\ell,u]$로 분할하고, 반복마다 고정된 KKT 행렬의 선형계 하나, 상자 투영, 승수 갱신만 수행한다. $\rho$가 바뀌지 않는 동안 분해는 한 번이므로 반복당 비용이 작고, 대신 반복 수가 $\rho$와 문제의 스케일에 크게 좌우된다.

소스(auxil.c)의 우변과 갱신식은 다음 반복에 해당하며, 이 장의 구현이 OSQP와 반복 단위로 일치함을 1.13.5 ②에서 확인했다. 첫 식은 KKT 선형계이며 $R=\operatorname{diag}(\rho_i)$이다.

$$\begin{bmatrix}P+\sigma I&A^\top\\A&-R^{-1}\end{bmatrix}\begin{bmatrix}\tilde x\\\nu\end{bmatrix}=\begin{bmatrix}\sigma x^k-q\\z^k-R^{-1}y^k\end{bmatrix},\qquad \tilde z=z^k+R^{-1}(\nu-y^k)$$

$$x^{k+1}=\alpha\tilde x+(1-\alpha)x^k,\quad z^{k+1}=\Pi_{[\ell,u]}\bigl(w+R^{-1}y^k\bigr),\quad y^{k+1}=y^k+R\,(w-z^{k+1}),\quad w=\alpha\tilde z+(1-\alpha)z^k$$

$\nu$를 소거하면 $(P+\sigma I+A^\top RA)\tilde x=\sigma x^k-q+A^\top(Rz^k-y^k)$, $\tilde z=A\tilde x$이다. OSQP는 소스의 기본 설정에서 직접 분해 해법기를 사용하고, 이 장의 코드는 위 응축형을 숄레스키로 푼다. 소스에서 확인한 규칙은 다음과 같다.

• 행별 $\rho$: 등식 행(상한과 하한의 차가 $10^{-4}$ 미만)은 $10^{3}\rho$, 나머지는 $\rho$이다. 소스의 주석은 제약 종류 $-1$에 OSQP_RHO_MIN($10^{-6}$)을 대응시키며 경계가 모두 무한한 행으로 판단한다 (확인 필요: 분류 함수).
• $\rho$ 적응: $\rho^{+}=\rho\sqrt{\dfrac{r_{\rm prim}/\max(\lVert Ax\rVert_\infty,\lVert z\rVert_\infty)}{r_{\rm dual}/\max(\lVert Px\rVert_\infty,\lVert A^\top y\rVert_\infty,\lVert q\rVert_\infty)}}$이며, 새 값이 현재의 5배를 넘거나 1/5 미만일 때만 채택하고 KKT 행렬을 다시 분해한다. 기본은 50회 반복마다 검사한다.
• 정지 기준: 25회마다 $r_{\rm prim}<\varepsilon_{\rm abs}+\varepsilon_{\rm rel}\max(\lVert Ax\rVert_\infty,\lVert z\rVert_\infty)$, $r_{\rm dual}<\varepsilon_{\rm abs}+\varepsilon_{\rm rel}\max(\lVert Px\rVert_\infty,\lVert A^\top y\rVert_\infty,\lVert q\rVert_\infty)$를 검사하고, 기본 설정에서는 쌍대성 간격도 같은 형식으로 검사한다.
• 불가능성 검출: $\delta y=y^{k+1}-y^k$에 대해 $\lVert A^\top\delta y\rVert_\infty<\varepsilon_{\rm prim\_inf}\lVert\delta y\rVert_\infty$이고 $S(\delta y)<0$이면 원 불가능으로 판정한다(Banjac 외 2019의 기준 (확인 필요: 논문의 정리)).
• 스케일링과 폴리싱: 기본은 평형화(Ruiz) 10회이다. 폴리싱은 $z_j-\ell_j<-y_j$이면 하한 활성, $u_j-z_j<y_j$이면 상한 활성으로 추정한 집합을 등식으로 둔 축소 KKT를 풀어 해를 정제하며 기본은 꺼져 있다.

위 표는 osqp_api_constants.h(master)의 값이며 pip 1.1.3 실행에서 같은 값을 읽었다.

1.13.3.6 근접 증강 라그랑주 — ProxQP

결론: ProxQP는 QP를 근접 곱셈자법의 외부 반복으로 풀고, 각 외부 반복의 부문제를 준매끄러운 뉴턴법으로 푼다. 뉴턴 단계는 작업 집합에 해당하는 행만 포함한 선형계이므로 활성 집합법과 닮았고, 근접 항 덕분에 행렬이 항상 정칙이다.

논문(식 8~10, 22)의 구조는 다음과 같다. 등식 $Ax=b$와 부등식 $l\le Cx\le u$에 대해 근접 항 $\rho,\mu_{eq},\mu_{in}$을 두고 $x^{k+1}\approx\arg\min\Phi_\rho^k(x)$, $y^{k+1}=y^k+(Ax^{k+1}-b)/\mu_{eq}$, $z^{k+1}=[z^k+(Cx^{k+1}-u)/\mu_{in}]_+$를 갱신한다. 부문제의 뉴턴 단계는 $\begin{bmatrix}H+\rho I&A^\top&C_{\mathcal I}^\top\\A&-\mu_{eq}I&0\\C_{\mathcal I}&0&-\mu_{in}I\end{bmatrix}$를 푼다. $\mathcal I$는 현재 활성인 부등식 행이다. 벌점 인자는 BCL 전략으로 갱신하고, 초기값은 등식만 둔 QP의 해(식 36)이며, 사전조건은 Ruiz 평형화이다. 기본 정지 조건은 원 잔차와 쌍대 잔차이고 쌍대성 간격은 check_duality_gap을 켤 때만 판정한다(문서 3-ProxQP_solve.md). 초기 추측 옵션은 NO_INITIAL_GUESS·EQUALITY_CONSTRAINED_INITIAL_GUESS·WARM_START_WITH_PREVIOUS_RESULT·WARM_START·COLD_START_WITH_PREVIOUS_RESULT의 다섯 가지이며, 근접 인자와 분해를 이어받는 것은 WARM_START_WITH_PREVIOUS_RESULT뿐이고 나머지는 기본값으로 초기화한다(문서의 설정 표). 논문은 비교 실험에서 자신의 방법이 작은 문제와 중간 크기 문제에서 효율·정확도·강건성이 좋다고 주장한다 (공개 주장치).

1.13.3.7 로봇의 QP 정식화

결론: 로봇의 세 대표 QP는 행렬 $P$·$A$가 시간에 따라 거의 변하지 않고 $q$·$\ell$·$u$만 바뀌는 구조이다. 이 구조 때문에 분해 재사용과 열 시작이 가능하다.

관절 토크는 $\tau=J_c^\top f$로 얻는다(12.7). 마찰 피라미드 $\lvert f_x\rvert\le\mu f_z$, $\lvert f_y\rvert\le\mu f_z$는 마찰 원뿔 $\sqrt{f_x^2+f_y^2}\le\mu f_z$를 외접하는 근사이며 $\mu/\sqrt2$로 줄인 피라미드가 내접 근사이다(마찰 원뿔은 13.2, 3.12). 응축형 선형 MPC는 예측 $X=\Phi x_0+\Gamma U$를 대입하여 $P=2(\Gamma^\top\bar Q\Gamma+\bar R)$, $q=2\Gamma^\top\bar Q(\Phi x_0-x_{\rm ref})$를 얻고, 상태 제약은 $\ell-\Phi x_0\le\Gamma U\le u-\Phi x_0$로 적힌다. 초기 상태 $x_0$와 목표가 바뀌어도 $P$와 $A=[I;\Gamma]$는 같다.

1.13.4 알고리즘

이 장의 구현 세 가지는 모두 robotics/fig-src/1.13-qp.py에 있다(dual_as·admm_step·ipm). 의사코드는 구현과 같은 순서이다.

쌍대 활성 집합법(Goldfarb–Idnani 형, $a_i^\top x\ge b_i$)

W ← 등식 제약 ∪ W0 ;  (x, μ) ← W를 등식으로 둔 QP의 해
μ_j < 0인 비등식 j를 하나씩 제거하며 (x, μ) 재계산        (열 시작 정리)
반복:
    p ← a_i^T x − b_i가 가장 작은 제약 ;  위반이 없으면 종료
    반복:
        [P A_W^T; A_W 0] [z; r] = [a_p; 0]               (원 방향 z, 승수 방향 r)
        t1 ← min{ μ_j / r_j : r_j > 0, j는 W의 비등식 }    (쌍대 단계 한계)
        t2 ← (b_p − a_p^T x) / (a_p^T z)                  (원 단계, a_p^T z ≤ 0이면 ∞)
        t = min(t1, t2) = ∞이면 불가능
        x ← x + t z ;  μ_W ← μ_W − t r ;  μ_p ← μ_p + t
        t = t2이면 p를 W에 더하고 바깥 반복으로 ;  아니면 t1의 j를 W에서 뺀다

OSQP 방식 ADMM

F ← chol(P + σI + A^T diag(ρ) A)                         (ρ가 바뀔 때만)
반복 k = 1, 2, …
    r ← σx − q + A^T(ρ∘z − y) ;  x̃ ← F⁻¹ r ;  z̃ ← A x̃
    x ← αx̃ + (1−α)x ;  w ← αz̃ + (1−α)z
    z⁺ ← clip(w + y/ρ, ℓ, u) ;  y ← y + ρ∘(w − z⁺) ;  z ← z⁺
    k가 25의 배수 : r_prim, r_dual, 쌍대성 간격이 모두 허용오차 이내이면 종료
    k가 50의 배수 : ρ ← ρ √(…), 5배 초과·1/5 미만이면 채택하고 F를 다시 분해

원-쌍대 내점법(Mehrotra)

(x, s, λ, ν) 초기화 ;  μ ← s^T λ / m
반복:
    r_d, r_e, r_p 계산 ;  모두 허용오차 이하이고 μ도 이하이면 종료
    D ← diag(λ/s) ;  [P + G^T D G, E^T ; E, 0] 분해
    예측 단계(σ=0) → α_aff, μ_aff ;  σ ← (μ_aff / μ)^3
    수정 단계(2차 항 Δs_aff∘Δλ_aff 포함, 중심화 σμ)
    α ← min(1, 0.99 · 경계까지의 최대 단계) ;  (x, s, λ, ν) ← (x, s, λ, ν) + αΔ

1.13.5 풀이 예제 — 세 해법의 구현과 라이브러리 비교

세 해법을 순수 파이썬으로 구현해 KKT 잔차와 반복 수를 측정하고, 같은 입력을 OSQP·ProxQP·qpOASES 등으로 풀어 대조한다. 결론부터 적으면 다음과 같다. ① 자체 쌍대 활성 집합법은 문제 B·C에서 KKT 잔차 $10^{-12}$ 이하의 해를 유한 반복에 준다. ② 자체 ADMM은 스케일링·폴리싱을 끈 OSQP와 반복 수가 같고 $x$가 $10^{-11}$ 이내로 같다. ③ $\rho$는 반복 수를 크게 바꾸며 문제 C에서 $\rho=0.1$은 10,000회에도 $10^{-6}$에 못 미친다. ④ 열 시작은 활성 집합법에서 평균 19.5회를 4.1회로 줄이지만 ADMM에서는 136.0회를 110.5회로 줄인다. ⑤ 반환 상태 '해결됨'은 최적을 보증하지 않는다. ProxQP는 기본 설정에서 무작위 QP 200개 가운데 19개에 대해 목적값 상대 오차 $10^{-3}$을 넘겼고(최대 1.2), OSQP는 기본 허용오차 $10^{-3}$ 수준의 제약 위반을 52개에서 보였다.

문제 A는 관절각 $(0.3,0.9,-0.6)$ rad의 3R 팔(1.1.5)에서 말단 속도 $(0,0.25)$ m/s를 요구하고 $\lambda=10^{-3}$, $\lvert\dot q_i\rvert\le0.5$ rad/s이다. 문제 B는 네 발이 $(\pm0.25,\pm0.15,-0.4)$ m에 있는 12 kg 몸체의 접촉력 배분이고 가속도 $(4,1,0)$ m/s², $\mu=0.5$, $0\le f_z\le60$ N이다. 문제 C는 이중 적분기($\Delta t=0.1$ s, $N=20$)로 $\lvert a\rvert\le2$ m/s², $\lvert v\rvert\le1$ m/s, $p\le1.5$ m이며 목표 위치 2.0 m, 위치·속도·입력 가중 $(10,1,0.1)$이다. 문제 T는 도해용이며 $\min\tfrac12(x-c)^\top H(x-c)$, $c=(3.2,2.9)$, $H=\left[\begin{smallmatrix}2&0.8\\0.8&1.5\end{smallmatrix}\right]$에 $x_1\le2.2$, $x_2\le2.4$, $x_1+x_2\le3.6$, $-x_1+2x_2\le2.6$, $x\ge0$이다. $\kappa(P)$는 1.1의 SVD로 계산했다.

ADMM의 반복 1회와 분해의 코드는 다음과 같다. 숄레스키와 삼각 풀이는 1.2의 구현이며, 정지 판정과 $\rho$ 갱신을 포함한 구동 함수 solve_admm과 dual_as·ipm은 스크립트에 있다.

def mv(M, v):
    return [sum(a * b for a, b in zip(r, v)) for r in M]

def mtv(M, v):
    return [sum(M[i][j] * v[i] for i in range(len(M)))
            for j in range(len(M[0]))]

def admm_factor(P, A, eqr, rho, sigma=1e-6):
    """K = P + sI + A^T R A의 숄레스키: (rv, L, L^T)"""
    n, m = len(P), len(A)
    rv = [1e3 * rho if e else rho for e in eqr]
    K = [[P[i][j] + (sigma if i == j else 0.0)
          + sum(rv[k] * A[k][i] * A[k][j] for k in range(m))
          for j in range(n)] for i in range(n)]
    L = chol(K)
    return rv, L, [list(c) for c in zip(*L)]

def admm_step(F, P, q, A, l, u, s, alpha=1.6, sigma=1e-6):
    """OSQP ADMM 1회: s=(x, z, y) -> 다음 (x, z, y)"""
    rv, L, Lt = F
    x, z, y = s
    m, n = len(l), len(q)
    r = mtv(A, [rv[i] * z[i] - y[i] for i in range(m)])
    r = [sigma * x[i] - q[i] + r[i] for i in range(n)]
    xt = tri_solve(Lt, tri_solve(L, r, True), False)
    zt = mv(A, xt)
    xn = [alpha * a + (1 - alpha) * b for a, b in zip(xt, x)]
    zr = [alpha * a + (1 - alpha) * b for a, b in zip(zt, z)]
    zn = [min(max(zr[i] + y[i] / rv[i], l[i]), u[i])
          for i in range(m)]
    yn = [y[i] + rv[i] * (zr[i] - zn[i]) for i in range(m)]
    return xn, zn, yn

① 해의 정확도와 반복 수. 반복 수는 $\max(r_{\rm prim},r_{\rm dual})\le10^{-6}$(ADMM, 매 반복 검사), 내점법은 모든 잔차와 $\mu$가 $10^{-6}$ 이하가 되는 반복이다.

쌍대 활성 집합법의 반복은 B에서 제약 3개 추가, C에서 25개 추가·10개 제거이며 최적 활성 집합은 C에서 입력 상한 5개·속도 상한 9개·위치 상한 1개이다. 내점법의 상보 잔차는 문제 C에서 $3.0\times10^{-9}$로 쌍대 잔차보다 크고 $x$ 오차는 $3.6\times10^{-7}$이다. 문제 A의 해는 관절 2·3의 속도 한계가 활성이고 말단 속도가 $(-0.067,0.193)$ m/s로 요구 $(0,0.25)$ m/s를 충족하지 못한다. 문제 B의 해는 법선력이 $(3.3,17.2,39.6,57.7)$ N(합 117.72 N)이고 마찰 제약 3개(첫째 발의 두 면, 둘째 발의 한 면)가 활성이다. 그림 1은 문제 T에서 세 계열이 같은 꼭짓점 $(2.2,1.4)$로 가는 경로이다.

② OSQP와의 반복 일치. 자체 ADMM과 OSQP 1.1.3(스케일링·폴리싱 끔)을 같은 입력과 설정으로 비교했다. 반복 수가 같고 $x$가 $10^{-11}$ 이내이다(문제 C의 4,000회 포함, $y$ 차이는 최대 $3.6\times10^{-11}$).

③ $\rho$ 민감도와 수렴 곡선. 반복 수는 $\rho$에 대해 U자형이며 최적 $\rho$는 A에서 0.01, C에서 10으로 세 자릿수 다르다. 적응 규칙은 최적 고정값의 1.1~1.3배 반복으로 수렴한다(A 33 대 28, B 67 대 52, C 837 대 646, T 113 대 104). 문제 C에서 $\max(r_{\rm prim},r_{\rm dual})$이 $10^{-6}$·$10^{-9}$·$10^{-12}$ 아래로 내려가는 반복은 $\rho=10$에서 646·1,040·1,364, 적응에서 837·1,332·1,792이고 $\rho=0.1$은 3,000회에 $1.4\times10^{-2}$이다. 내점법은 11회에 $10^{-10}$ 수준이다.

④ 열 시작. 문제 C의 폐루프(50단계, 단계 25에서 목표가 2.0 m에서 0 m로 바뀐다)에서 쌍대 활성 집합법의 해로 진행하며 얻은 QP 50개를 같은 순서로 풀었다. 냉 시작(cold start)은 해·승수를 0으로, 작업 집합을 등식 제약만으로 시작한 경우이며 OSQP는 warm_starting을 끈 경우, ProxQP는 NO_INITIAL_GUESS, qpOASES는 단계마다 새 함수 객체(init)이다. 활성 집합법의 열 시작은 24단계에서 반복이 0회이고 목표가 바뀐 단계 25에서 22회이다. 반복의 정의가 해법기마다 다르므로 표의 값은 같은 해법기의 냉·열 시작 비교에만 사용한다.

OSQP의 스케일링을 끄면 평균 반복이 56.0에서 136.0으로 늘어, 이 문제에서 스케일링이 열 시작보다 큰 효과를 가진다. ProxQP에서 근접 인자를 유지하는 WARM_START_WITH_PREVIOUS_RESULT는 16개 단계에서 원 불가능을 반환했고 그 가운데 한 단계는 10,006회에 이르렀다(근접 인자 $\mu_{in}$이 단계마다 $10^{-3}$에서 $10^{-8}$로 줄었다). 근접 인자를 초기화하는 COLD_START_WITH_PREVIOUS_RESULT는 반복을 줄이지 못했다. qpOASES는 CasADi의 인터페이스가 같은 함수 객체의 두 번째 호출부터 hotstart를 부르는 것을 소스에서 읽었고(called_once), 평균 nWSR이 46.7에서 4.8로 줄었다.

⑤ 라이브러리 비교. 같은 문제 B·C의 해를 활성 집합법의 기준해와 비교했다. 반복 칸은 해법기마다 정의가 다르다(활성 집합법: 추가·제거 수, ADMM: 반복, 내점법: 뉴턴 단계, qpOASES: nWSR, ProxQP: 내부 뉴턴 반복). 소요 시간은 호출 오버헤드·기본 허용오차·반복의 정의가 달라 같은 조건이 아니므로 비교하지 않았다.

quadprog(Goldfarb–Idnani 구현)는 문제 C에서 21회 추가·5회 제거로 자체 구현의 25·10과 다르다. 위반 제약의 선택 규칙의 차이로 추정하며 확인하지 않았다 (추정). 기본 설정의 정확도 차이는 상당 부분 기본 허용오차에서 온다. OSQP 기본($\varepsilon=10^{-3}$)의 문제 C 해는 $x$ 오차가 $4.9\times10^{-2}$(가속도 한계 2의 2.4%)이고 $\varepsilon=10^{-6}$에서 $2.1\times10^{-5}$이다. CVXPY는 OSQP를 선택하며 소스(osqp_qpif.py)에서 $\varepsilon_{\rm abs}=\varepsilon_{\rm rel}=10^{-5}$, 최대 반복 10,000, 폴리싱 켜짐으로 기본값을 바꾼다.

⑥ 상태 코드의 신뢰도. 무작위 볼록 QP 200개(시드 31, 변수 3~9개, 부등식 $n$~$2n+1$개, 모두 가능)와 위 QP 수열 50개를 각 해법기로 독립으로 풀었다. 첫 표의 '초과'는 목적값 상대 오차 $(f(x)-f^\ast)/\max(1,\lvert f^\ast\rvert)$ 또는 제약 위반이 $10^{-3}$을 넘은 개수이고, '상태 불량'은 해법기가 해결 상태를 반환하지 않은 개수이다. 둘째 표는 QP 수열에서 해결 상태를 반환하지 않은 횟수와 단계, 나머지 단계의 최대 $x$ 오차이다.

OSQP $\varepsilon=10^{-6}$은 단계 21~24에서 최대 반복 4,000회에 도달했다(maximum iterations reached). ProxQP의 기본 정지 조건에는 쌍대성 간격이 없다. 소스는 원·쌍대 잔차가 허용오차 이내이면 종료하며, 문서는 이 조건만으로 KKT 조건이 모두 만족됨을 보증하지 않고 check_duality_gap이 기본적으로 꺼져 있다고 밝히며 다른 해법기도 일반적으로 이 기준을 검사하지 않는다고 설명한다(OSQP 1.0 소스의 기본은 검사한다). 무작위 QP의 최악 사례(변수 8개·부등식 9개)에서 기본 설정은 3회 반복 뒤 '해결됨'으로 종료했고 반환 목적값은 0.324로 최적값 $-1.342$와 달랐으며 보고된 쌍대성 간격은 0.43이었다. 기본 설정의 반복은 평균 7.7회이고 17개 QP는 2회 이하에서 종료했다. 위 QP 수열에서는 기본 설정이 가능한 QP 7개(단계 18~24)를 PROXQP_PRIMAL_INFEASIBLE로 판정했다. 이 QP들의 체비쇼프 여유는 0.15~1.0이므로 가능 영역은 내부를 가지며, eps_primal_inf를 $10^{-9}$로 낮추거나 NO_INITIAL_GUESS로 시작하면 7개 모두 해결로 반환했다. 원인은 규명하지 못했다. 소스에서 원 불가능성 판정은 내부 뉴턴 반복마다 뉴턴 방향 $(\delta y,\delta z)$로 수행되며(frequence_infeasibility_check=1) 판정식은 읽지 않았다.

⑦ 불가능한 문제. 문제 B에서 가속도를 $(6,1,0)$ m/s²로 높이면 마찰 피라미드가 요구 수평력 72 N을 허용하는 한계 약 59 N(총 법선력 117.7 N × 0.5)을 넘어 불가능하다.

OSQP는 75회(스케일링 끔 100회)에 원 불가능을 판정하고 반환한 $x$는 의미가 없다. quadprog는 예외를 던지고 qpOASES는 CasADi의 error_on_fail을 끄면 상태 문자열을 반환하며 켜 둔 기본에서는 예외(RuntimeError)를 던졌다. 호출부는 불가능 판정을 정상 흐름으로 처리해야 한다(1.13.9).

한계. 수치는 합성 입력(시드 고정)이며 문제 4종과 무작위 QP로 측정했다. 자체 구현은 밀집 선형대수이고 스케일링이 없어 성능 비교의 대상이 아니다. 라이브러리는 pip 바이너리를 기본 설정에 가까운 상태로 실행했으며 해법기마다 최적 설정을 찾지 않았다. qpOASES는 CasADi 3.8.2에 번들된 것이며 저작권 배너의 연도가 2007–2015이고 버전은 식별하지 못했다 (확인 필요).

1.13.6 설계 트레이드오프 — 해법기의 선택

1.13.6.1 세 계열의 비교

견해: 이 장의 측정 범위(변수 2~20개의 밀집 QP, 같은 구조의 QP 수열)에서 제어 주기마다 비슷한 QP를 푸는 MPC·전신 제어에는 열 시작이 효과적인 해법을 우선 평가하고, 최악 실행 시간이 중요하면 반복 수 상한을 두고 상한 도달 시의 대체 동작을 설계하는 것이 해법기 선택보다 중요하다고 본다. 근거는 활성 집합법의 열 시작 반복이 첫 QP를 제외하면 최대 22회(목표 변경 단계)였고, ADMM의 평균 반복이 스케일링 유무로 56.0에서 136.0으로 달라졌으며, OSQP·ProxQP는 설정에 따라 상태가 불량해졌다는 측정이다. 대형 희소 문제에서는 다르게 나올 수 있다.

1.13.6.2 대표 해법기와 현황

기준 시점 2026-11-03. 제원은 GitHub API의 공개 표기이며 푸시 날짜는 유지의 대리 지표일 뿐이다(0.16). 상용 해법기는 이 장에서 다루지 않았다.

qpOASES는 릴리스 없이 태그만 있으며(최신 태그 releases/3.2.2) 읽은 매뉴얼은 3.2판(2017-04)이다. OSQP 소스에는 임베디드 모드(OSQP_EMBEDDED_MODE)에서 $\rho$ 갱신 코드를 제외하는 조건부 컴파일이 있고 파이썬 패키지에는 코드 생성(codegen) 메서드가 있다(1.15). 해법기의 기본 정지 조건은 다음과 같다.

1.13.7 응용 축 비교

이 표는 이 책의 정리이며 장별 요구는 다르다.

1.13.8 관련 규격·프로토콜

QP 해법기의 입력·출력 형식을 정하는 국제 표준은 없다. 해법기마다 입력 형태($\ell\le Ax\le u$, $Ax=b$와 $l\le Cx\le u$ 등)와 부호 규약이 다르므로 사용하는 쪽의 문서를 읽는 것이 상호운용의 기준이다.

1.13.9 관련 규제·인증 요구와의 연결

QP 해법기 자체에 적용되는 규제는 없다. 해법기가 속도·힘 한계의 감시나 안전 필터(21.4)처럼 안전 기능의 일부이면 소프트웨어 안전 요구(26.20)와 기능안전 표준(26.5, 차량은 27.15)이 그 입력·출력과 실행 시간에 적용된다. 이 장의 측정에서 도출되는 검증 항목은 세 가지이다. 상태 코드와 별개로 KKT 잔차를 계산해 해를 확인할 것, 불가능 판정과 최대 반복 도달을 정상 흐름으로 처리하고 대체 동작을 정의할 것, 허용오차·초기 추측·스케일링 설정을 형상 관리 대상으로 둘 것이다.

1.13.10 한계·유의점·미해결

• 자체 구현의 범위. 밀집 선형대수, 스케일링 없음, 퇴화 처리 없음, 쌍대 활성 집합법은 단계마다 KKT 행렬을 새로 푼다(분해 갱신 없음). 쌍대 활성 집합법은 $P\succ0$을 가정한다. 희소·대형·GPU 구현의 거동은 측정하지 않았다.
• 라이브러리 실행의 범위. pip 바이너리의 기본에 가까운 설정이며 해법기별 튜닝은 하지 않았다. Clarabel·CVXPY·SciPy는 비교 목적의 참고 행이다. 시간 측정은 하지 않았다.
• 읽지 못한 것. OSQP 논문과 ADMM 수렴 이론, Goldfarb–Idnani·Mehrotra의 원 논문, qpOASES 2014년 논문은 서지만 확인했다. 소스·매뉴얼에서 읽은 범위를 넘는 알고리즘 성질은 (확인 필요)로 표시했다.
• ProxQP의 관측. 가능한 QP에 대한 원 불가능 판정과 기본 정지 조건의 최적성 손실은 이 장의 문제군과 설정에서 재현되었다. 원인은 규명하지 못했고 다른 문제군과 이후 버전에서의 거동은 확인하지 않았다.
• 불가능·퇴화. 퇴화 QP에서 활성 집합법의 순환 방지와 ADMM의 불가능 인증은 측정하지 않았다(불가능 사례 1건만 확인).
• 비볼록·정수. 비볼록 QP와 혼합정수 QP는 1.14, 모델링·코드 생성은 1.15, 다단계 계층 QP는 12.7이 다룬다.

1.13.11 체크포인트·연구 질문

체크포인트:

1. $\ell\le Ax\le u$ 형태로 등식·경계·일방 제약을 적고 KKT 조건과 쌍대성 간격을 계산할 수 있는가.
2. 작업 집합의 등식 QP에서 승수의 부호로 제약을 제거하는 이유를 설명할 수 있는가.
3. 쌍대 활성 집합법의 단계 길이 $t_1$, $t_2$의 의미를 설명할 수 있는가.
4. 내점법의 중심 경로와 반복마다 분해를 다시 하는 이유를 설명할 수 있는가.
5. OSQP의 반복식에서 $\rho$·$\sigma$·$\alpha$의 역할과 $\rho$ 적응·스케일링이 반복 수에 주는 영향을 설명할 수 있는가.
6. 열 시작이 활성 집합법·ADMM·내점법에서 각각 얼마나 효과적인지 근거와 함께 말할 수 있는가.
7. 해법기가 반환한 '해결됨'을 신뢰하기 전에 확인할 세 가지(KKT 잔차·쌍대성 간격·설정)를 말할 수 있는가.

연구 질문(47.12에 수록한다):

• 로봇 제어용 QP 수열에서 반복 수·지연의 분포(최악 경우 포함)를 해법기·설정별로 측정한 공개 데이터는 어느 정도이며, 반복 수 상한을 둔 실시간 보증은 어떻게 검증되는가.
• 기본 정지 조건(쌍대성 간격·상보성 미판정, 낮은 허용오차)이 로봇 제어 QP의 최적성 손실로 이어지는 빈도와 크기는 얼마인가.
• 가능한 QP에 대한 원 불가능 판정의 오탐을 줄이는 초기 추측·판정 규칙은 무엇인가.
• 스케일링과 $\rho$ 선택을 문제 구조(MPC·접촉력 배분)에 맞춰 자동화하면 반복 수의 분산이 얼마나 줄어드는가.

1차 출처

• S. Boyd, L. Vandenberghe, 『Convex Optimization』, Cambridge University Press, 2004(doi:10.1017/CBO9780511804441) (확인 필요: 장 번호).
• J. Nocedal, S. J. Wright, 『Numerical Optimization』, 2판, Springer, 2006(doi:10.1007/978-0-387-40065-5) — 16장 'Quadratic Programming'(Crossref 장 항목으로 제목 확인)은 본문 미확인.
• D. Goldfarb, A. Idnani, 「A numerically stable dual method for solving strictly convex quadratic programs」, Mathematical Programming 27(1), 1–33, 1983(doi:10.1007/BF02591962) — 서지만 확인.
• S. Mehrotra, 「On the implementation of a primal-dual interior point method」, SIAM Journal on Optimization 2(4), 575–601, 1992(doi:10.1137/0802028) — 서지만 확인.
• B. Stellato, G. Banjac, P. Goulart, A. Bemporad, S. Boyd, 「OSQP: an operator splitting solver for quadratic programs」, Mathematical Programming Computation 12(4), 637–672, 2020(doi:10.1007/s12532-020-00179-2; arXiv:1711.08013) — 서지만 확인(arXiv API 대조), 소스는 읽음.
• G. Banjac, P. Goulart, B. Stellato, S. Boyd, 「Infeasibility detection in the alternating direction method of multipliers for convex optimization」, Journal of Optimization Theory and Applications 183(2), 490–519, 2019(doi:10.1007/s10957-019-01575-y) — 서지만 확인.
• S. Boyd, N. Parikh, E. Chu, B. Peleato, J. Eckstein, 「Distributed optimization and statistical learning via the alternating direction method of multipliers」, Foundations and Trends in Machine Learning 3(1), 1–122, 2011(doi:10.1561/2200000016) — 서지만 확인.
• H. J. Ferreau, C. Kirches, A. Potschka, H. G. Bock, M. Diehl, 「qpOASES: a parametric active-set algorithm for quadratic programming」, Mathematical Programming Computation 6(4), 327–363, 2014(doi:10.1007/s12532-014-0071-1) — 서지만 확인, 사용자 매뉴얼은 읽음.
• H. J. Ferreau, H. G. Bock, M. Diehl, 「An online active set strategy to overcome the limitations of explicit MPC」, International Journal of Robust and Nonlinear Control 18(8), 816–830(매뉴얼 표기 2008, Crossref 발행 연도 2007; doi:10.1002/rnc.1251) — 서지만 확인.
• A. Bambade, S. El-Kazdadi, A. Taylor, J. Carpentier, 「PROX-QP: Yet another quadratic programming solver for robotics and beyond」, RSS 2022(HAL hal-03683733) — HAL 판을 읽음.
• R. T. Rockafellar, 「Augmented Lagrangians and applications of the proximal point algorithm in convex programming」, Mathematics of Operations Research 1(2), 97–116, 1976(doi:10.1287/moor.1.2.97) — 서지만 확인.
• F. Borrelli, A. Bemporad, M. Morari, 『Predictive Control for Linear and Hybrid Systems』, Cambridge University Press, 2017(doi:10.1017/9781139061759) — 서지만 확인.
• I. Maros, C. Mészáros, 「A repository of convex quadratic programming problems」, Optimization Methods and Software 11(1–4), 671–681, 1999(doi:10.1080/10556789908805768) — 서지만 확인.
• P. J. Goulart, Y. Chen, 「Clarabel: An interior-point solver for conic programs with quadratic objectives」, arXiv:2405.12762(Mathematical Programming Computation 등록 2026, doi:10.1007/s12532-026-00320-7 (확인 필요: 권·쪽)) — 서지만 확인.
• D. Escande, N. Mansard, P.-B. Wieber, 「Hierarchical quadratic programming: Fast online humanoid-robot motion generation」, International Journal of Robotics Research 33(7), 1006–1028, 2014(doi:10.1177/0278364914521306) — 서지만 확인(계층 QP는 12.7).
• IEEE Std 754-2019, 『IEEE Standard for Floating-Point Arithmetic』(doi:10.1109/IEEESTD.2019.8766229) — 서지만 확인.
• 소스·문서(조회·실행 2026-11-03): osqp/osqp(osqp_api_constants.h·osqp_api.c·auxil.c·polish.c, master)와 osqp 1.1.3 실행, qpOASES 3.2 사용자 매뉴얼·README·Options.cpp(coin-or/qpOASES), Simple-Robotics/proxsuite(settings.hpp·solver.hpp·doc/2-ProxQP_api.md·doc/3-ProxQP_solve.md, devel)와 proxsuite 0.8.2 실행, casadi/casadi의 qpoases_interface.cpp와 casadi 3.8.2 실행, cvxpy 1.9.3의 osqp_qpif.py, quadprog 0.1.14(PyPI 메타데이터), clarabel 0.11.1·scipy 1.18.1 실행. 실행 스크립트는 robotics/fig-src/1.13-libcheck.py, 출력은 robotics/fig-src/data/1.13-libcheck-output.txt와 1.13-libcheck-repos-output.txt에 있다.

다음 장. 1.14는 최적화 III으로, 비선형 계획법의 SQP·내점법(IPOPT·SNOPT)과 혼합정수 계획법(Gurobi)을 다룬다. 이 장의 QP는 SQP의 부문제로 반복해서 풀리며, 활성 집합과 열 시작의 성질이 비선형 계획법의 수렴과 비용을 결정한다.