以下是一種循序推理的方式，來幫助你從基礎(chǔ)概念出發(fā)，理解 內(nèi)點(diǎn)法（Interior-Point Method, IPM） 是什么、為什么要用它，以及它是如何工作的。

---

1. 問題起點(diǎn)：帶不等式約束的優(yōu)化

假設(shè)你有一個(gè)帶不等式約束的優(yōu)化問題：
$$\begin{array}{rl}& \underset{x}{\min }f(x)\\ & \text{subject?to?}{g}_i(x)\le 0,i=1,\dots ,m.\end{array}$$

• $f(x)$：目標(biāo)函數(shù)
• $g_i(x)\le 0$：不等式約束（例如，控制輸入不能超范圍）

我們想要找到一個(gè)滿足所有約束且使目標(biāo)函數(shù)最小的點(diǎn) $x^{?}$。

---

2. 直接在邊界上“走”的困難

一種思路是：最優(yōu)點(diǎn)往往在可行域邊界上（很多經(jīng)典優(yōu)化問題里，最優(yōu)解會(huì)出現(xiàn)在約束生效的邊界）。

• 但是，如果我們只在邊界上“走”，常常要先找到并跟隨所有活躍約束的交線，這在高維情況下非常復(fù)雜。
• 并且，如果問題是非線性的，直接處理“在邊界上走”會(huì)更難，容易出現(xiàn)數(shù)值不穩(wěn)定或無法收斂的問題。

---

3. 想法：能不能“在里面”搜索，最后再逼近邊界？

為了避免“一上來就卡在邊界”，有人提出：

• 先在可行域內(nèi)部（所有約束 $g_i(x)<0$ 都“寬?！?#xff09;附近做搜索，那里沒有太多束縛，比較容易用連續(xù)的梯度法去迭代尋優(yōu)；
• 當(dāng)我們逐漸找到一個(gè)更優(yōu)解，需要靠近邊界時(shí)，再“慢慢地”逼近它。

這樣做，就不需要每一步都嚴(yán)格跟在邊界上，也能保證可行性（不跑到約束外面去）。

---

4. 內(nèi)點(diǎn)法的關(guān)鍵：障礙函數(shù)（Barrier Function）

為在域內(nèi)部搜索，需要一種數(shù)學(xué)手段，讓解“自動(dòng)”不跑出界。這時(shí)就引入了障礙函數(shù)。

4.1 形式：對約束加一個(gè)“負(fù)對數(shù)”項(xiàng)

對每個(gè)不等式約束 $g_i(x)\le 0$，定義一個(gè)項(xiàng)：
$$?\log (?g_i(x)).$$

• 如果 $g_i(x)<0$，那么 $?g_i(x)>0$，對數(shù)是可以求的。
• 一旦 $g_i(x)$ 接近 0（也就是接近邊界），$?g_i(x)$ 趨向 0，這個(gè)對數(shù)會(huì)趨向 $?\infty$，但我們在目標(biāo)函數(shù)中是帶一個(gè)負(fù)號“-”，變成了一個(gè)正無窮大的懲罰。
• 這樣一來，就相當(dāng)于在可行域的邊緣設(shè)置了一道**“墻”，越靠近邊界，代價(jià)就越大，解被迫**留在可行域內(nèi)部。

4.2 組合成“障礙型”目標(biāo)函數(shù)

把原本的目標(biāo)函數(shù) $f(x)$ 和障礙項(xiàng)結(jié)合：
$$B(x,\mu )=f(x)?\mu \sum _{i=1}^m\log (?g_i(x)).$$

• $\mu >0$ 是一個(gè)系數(shù)，稱為“障礙參數(shù)（barrier parameter）”。
• 當(dāng) $\mu$ 較大時(shí)，障礙懲罰也大，解就離約束邊界更遠(yuǎn)；當(dāng)我們逐漸減小 $\mu$，系統(tǒng)會(huì)允許解更靠近邊界。

---

5. 迭代逼近最優(yōu)解

內(nèi)點(diǎn)法并不是一次就把問題解決，而是：

1. 從一個(gè)“較寬松”的問題開始（$\mu$ 比較大，邊界懲罰很強(qiáng)）。
2. 求解 ${\min }_{x}B(x,\mu )$，得到一個(gè)可行域內(nèi)部的解。
3. 降低 $\mu$ 的值，重新求解，這時(shí)可以允許解更加靠近（或到達(dá)）真正的最優(yōu)邊界。
4. 多次迭代后，$\mu \to 0$，解會(huì)逼近真實(shí)最優(yōu)解。

這一過程里，算法會(huì)用到類似于 牛頓法、KKT 條件 等工具來求解各輪的子問題。

---

6. 內(nèi)點(diǎn)法 vs. 其他方法

• 單純形法（Simplex）：在可行域的“頂點(diǎn)—邊界—面”上走，比較適合線性規(guī)劃；但在高維非線性問題上，效率可能較低。
• 內(nèi)點(diǎn)法：通過“障礙函數(shù)”把解留在域內(nèi)部，逐漸往邊界逼近，常在非線性規(guī)劃（NLP）中表現(xiàn)更好，也更適合大規(guī)模問題。

---

7. 在 MPC 中的應(yīng)用

MPC 求解器（如 IPOPT）正是利用內(nèi)點(diǎn)法來處理：

• 多步預(yù)測的系統(tǒng)動(dòng)力學(xué)約束
• 輸入/狀態(tài)上下限
• 非線性目標(biāo)函數(shù)

讓無人機(jī)等系統(tǒng)在高維空間里高效地搜索最優(yōu)控制輸入。

---

🔑 小結(jié)：內(nèi)點(diǎn)法的一句話總結(jié)

內(nèi)點(diǎn)法是一種在可行域“內(nèi)部”迭代搜索的求解策略，通過“障礙函數(shù)”阻擋解跑出可行域，并逐步放松障礙參數(shù)，最終逼近最優(yōu)解和約束邊界。

這就是內(nèi)點(diǎn)法的核心“推理過程”：與其從邊界開始，不如在“里層”走，讓數(shù)值算法更穩(wěn)定，再慢慢讓解貼近約束邊界，從而找到最優(yōu)。

下面我將針對這段代碼的邏輯和實(shí)現(xiàn)細(xì)節(jié)，結(jié)合“內(nèi)點(diǎn)法(障礙函數(shù)) + 固定步長梯度下降”這一思路，做一個(gè)比較細(xì)致的解析。

---

1. 整體思路與算法框架

目標(biāo)問題（從注釋和注解上看）是：

$$\begin{array}{rl}& \min f(x)=x_1+x_2+x_3\\ & \text{s.t.}{x}_1^2+x_2^2+x_3^2\le 1.\end{array}$$

• 可行域是單位球 { x ∈ R 3 : ∥ x ∥ ≤ 1 } \{x \in \mathbb{R}^3 : \|x\|\le 1 \} {x∈R3:∥x∥≤1}；

• 希望用內(nèi)點(diǎn)法來解決不等式約束 $x_1^2+x_2^2+x_3^2\le 1$。
  通常在內(nèi)點(diǎn)法中，會(huì)把不等式 (g(x) \le 0)（這里 $g(x)=x^2+y^2+z^2?1$ ）轉(zhuǎn)寫成形如 (h(x) > 0)，再將 (\log h(x)) 當(dāng)作障礙項(xiàng)（barrier）。本例里定義
  $$h(x)=1?(x_1^2+x_2^2+x_3^2),$$
  于是$h(x)>0$等價(jià)于 $x_1^2+x_2^2+x_3^2<1$。

• 內(nèi)點(diǎn)法的障礙目標(biāo)函數(shù)定義為
  $$B(x,\mu )=f(x)?\mu \ln (h(x)),$$
  其中 $\mu >0$ 是障礙系數(shù)(barrier parameter)。$\mu$ 越小，$?\mu \log h(x)$的懲罰作用越強(qiáng)，最終會(huì)將可行解“推”到最接近邊界的最優(yōu)位置上。

• 在固定一個(gè)$\mu$ 的情況下，常用牛頓法或梯度法去最小化 $B(x,\mu )$ 。在算法外層循環(huán)中，則逐步減小 $\mu$（比如乘個(gè)衰減因子），讓解逐步逼近約束邊界并最終得到較準(zhǔn)確的可行最優(yōu)解。

在這份代碼中：

1. 外層循環(huán) (共 max_outer 次)：

   • 每次先用當(dāng)前 $\mu$ 在球內(nèi)做若干步梯度下降，嘗試求解 $\min B(x,\mu )$；
   • 之后衰減 $\mu$ $\leftarrow$ $\mu ?$ $\text{(mu\_decay)}$，進(jìn)入下一輪；
   • 重復(fù)，直到$\mu$ 足夠小或者達(dá)到外層迭代上限。

2. 內(nèi)層循環(huán) (共 max_inner 次)：

   • 計(jì)算當(dāng)前 $B$ 的梯度 $?B(x,\mu )$；
   • 做一步固定步長的梯度下降：$x\leftarrow x?\alpha ?B(x,\mu )$；
   • 如果發(fā)現(xiàn)新點(diǎn) $x_{\text{new}}$不可行（或者離邊界過近），就縮小步長再試；
   • 如果兩次迭代 $\Vert x_{\text{new}}?x\Vert$ 非常小(< tol)，就視為收斂并 break。

這樣就得到一條漸進(jìn)接近最優(yōu)解的迭代軌跡，存放在 x_history 中。

---

2. 代碼主干解析

從最外層開始看起，核心部分在函數(shù)

function x_history = interior_point_3d_solve()% 參數(shù)mu_init = 1.0;      % 初始障礙參數(shù)mu_decay = 0.2;     % 每輪迭代后 mu 的衰減因子alpha    = 0.001;   % 固定梯度下降步長tol      = 1e-6;    % 收斂判據(jù)max_outer= 10;      % 外層循環(huán)次數(shù)max_inner= 50;      % 每次 mu 下最大內(nèi)層迭代次數(shù)% 初始化可行解(球內(nèi))x = [0;0;0];  mu = mu_init;x_history = [];for outer = 1:max_outerfor inner = 1:max_innerg = grad_B(x, mu);  x_new = x - alpha*g;% 若越過球邊界 => h(x_new) <= 0if h_3d(x_new) <= 1e-9x_new = x - 0.1*alpha*g;  % 縮小步長再試end% 收斂判斷if norm(x_new - x) < tolx = x_new;break;  % 跳出內(nèi)層循環(huán)endx = x_new;x_history = [x_history; x']; end% 減小 mumu = mu * mu_decay;if mu < 1e-12break;    % mu 已經(jīng)很小，不必再迭代endend% 加入最終點(diǎn)x_history = [x_history; x'];
end

1. mu_init 設(shè)為 1.0，之后每次外層循環(huán)會(huì) mu = mu * mu_decay，即乘以 0.2。這樣大概迭代幾次后就會(huì)讓 mu 變得很小。
2. alpha=0.001 是固定的梯度下降步長，相對比較小，所以我們用到 max_inner=50 步來讓它收斂到一個(gè)合適精度。
3. 收斂閾值 tol=1e-6：若相鄰兩步$\Vert x_{\text{new}}?x\Vert$ 小于這個(gè)值，就認(rèn)為內(nèi)層已經(jīng)收斂。
4. 每次更新 x 后都會(huì)把它記錄到 x_history 里，用于可視化迭代軌跡。

這里值得注意的是：如果單純用固定步長，可能碰到越過可行域邊界（即 (x^2+y2+z^2>1)）的風(fēng)險(xiǎn)。為此，代碼做了一個(gè)簡單檢查：

if h_3d(x_new) <= 1e-9x_new = x - 0.1*alpha*g;  % 步長縮小10倍
end

當(dāng)然，這只是一個(gè)很“簡單粗暴”的處理，工業(yè)級內(nèi)點(diǎn)法通常要做更精細(xì)的線搜索或牛頓校正，這里是為了示例演示。

---

3. 障礙函數(shù)與梯度的計(jì)算

3.1 h_3d(x)

function val = h_3d(x)val = 1.0 - (x(1)^2 + x(2)^2 + x(3)^2);
end

即$h(x)=1?r^2$，其中 $r^2=x_1^2+x_2^2+x_3^2$。球內(nèi)保證 $h(x)>0$。

3.2 B_3d(x, mu)

function val = B_3d(x, mu)val = f_3d(x) - mu*log(h_3d(x));
end

對應(yīng)障礙目標(biāo)$B(x,\mu )=f(x)?\mu \ln (h(x))$。

3.3 grad_B(x, mu) 的推導(dǎo)

從數(shù)學(xué)上看，如果
$$B(x,\mu )=f(x)?\mu \ln (h(x)),$$
則其梯度為
$$?B(x,\mu )=?f(x)?\mu ?\ln (h(x)).$$
其中
$$?\ln (h(x))=\frac{1}{h(x)}?h(x).$$
而 (h(x) = 1 - r^2)，(\nabla h(x) = -2x)。所以
$$?\ln (h(x))=\frac{?2x}{1?r^2}.$$
帶上負(fù)號“$?\mu$”一起，就得到對障礙項(xiàng)的貢獻(xiàn)為$+\frac{2\mu x}{1?r^2}$。

如果我們要最小化 $f(x)=x_1+x_2+x_3$，其梯度就是 $(1,1,1)$。
于是
$$?B(x,\mu )=(1,1,1)+\frac{2\mu }{1?r^2}(x_1,x_2,x_3).$$

在代碼中，可以看到：

function g = grad_B(x, mu)hx = h_3d(x);  % 1 - r^2dB_dx0 = 1.0 + (2.0 * mu * x(1) / hx);dB_dx1 = 1.0 + (2.0 * mu * x(2) / hx);dB_dx2 = 1.0 + (2.0 * mu * x(3) / hx);g = [dB_dx0; dB_dx1; dB_dx2];
end

正好對應(yīng)上面的公式：1.0 就是 $?f(x)$ 的那一部分，$(2.0?mu?x(i)/hx)$ 對應(yīng)障礙項(xiàng)梯度那一部分。

---

---

4. 運(yùn)行與結(jié)果

如果修正了 f_3d，那么在球內(nèi)最小化 (x+y+z) 的最優(yōu)解，理論上會(huì)落在球面上與 ((1,1,1)) 方向相反的地方，也就是球面朝著 ((-1,-1,-1)) 的方向。其最優(yōu)解應(yīng)當(dāng)是

$$\left(?\frac{1}{\sqrt{3}},?\frac{1}{\sqrt{3}},?\frac{1}{\sqrt{3}}\right)$$
因?yàn)樵诩s束 $x^2+y^2+z^2=1$ 上，(x+y+z) 的最小值就是 $?\sqrt{3}$。迭代跑起來后，你應(yīng)該能看到解最終趨近這個(gè)點(diǎn)，軌跡會(huì)從球心出發(fā)，一路在球內(nèi)前進(jìn)，并在迭代后期逐漸貼近球面。

如果把可視化部分加上（即下面這段示例）：

x_history = interior_point_3d_solve();figure('Color','w','Name','Interior-Point 3D Demo');
hold on; grid on; axis equal;% 繪制單位球面
[Xs, Ys, Zs] = sphere(50);
surf(Xs, Ys, Zs, ...'FaceAlpha',0.1, 'EdgeColor','none', 'FaceColor','c');xlabel('x_0'); ylabel('x_1'); zlabel('x_2');
title('Minimize x_0 + x_1 + x_2 subject to x_0^2 + x_1^2 + x_2^2 \le 1');% 畫迭代軌跡
x0_hist = x_history(:,1);
x1_hist = x_history(:,2);
x2_hist = x_history(:,3);
plot3(x0_hist, x1_hist, x2_hist, 'b-o','LineWidth',1.5);% 最后一點(diǎn)標(biāo)紅
plot3(x0_hist(end), x1_hist(end), x2_hist(end), ...'ro', 'MarkerSize',8, 'MarkerFaceColor','r');legend('Unit Sphere','Iter Process','Final Solution');
view(35,25);

就能看到一個(gè)球面和從原點(diǎn)出發(fā)，沿著負(fù)對角線方向慢慢收斂到球面那一點(diǎn)的軌跡。

---

5. 小結(jié)

1. 內(nèi)點(diǎn)法原理：通過把約束$x_1^2+x_2^2+x_3^2\le 1$ 轉(zhuǎn)化為對數(shù)障礙$?\mu \ln (1?r^2)$，并在外層迭代不斷減小 $\mu$。這能保證迭代點(diǎn)始終在球內(nèi) ($h(x)>0$)，同時(shí)在 $\mu \to 0$ 時(shí)逐漸收斂到可行域邊界上的最優(yōu)點(diǎn)。

2. 實(shí)現(xiàn)細(xì)節(jié)：

   • 用固定步長 + 簡單可行性校正(越界時(shí)縮步長)；
   • 用 max_outer 與 max_inner 控制多重循環(huán)；
   • 用 tol 判斷迭代收斂；
   • 在每次迭代都記錄 x 到 x_history 中，用于可視化。

3. 需要修正的地方：
   代碼中的 f_3d(x) 與注釋/梯度公式不一致，應(yīng)改回

   function val = f_3d(x)val = x(1) + x(2) + x(3);
   end
   

   這樣才與“最小化 (x+y+z)”的需求相吻合。

修正后運(yùn)行，即可得到一個(gè)從球心(原點(diǎn))出發(fā)，最終靠近 $(?\frac{1}{\sqrt{3}},?\frac{1}{\sqrt{3}},?\frac{1}{\sqrt{3}})$ 的迭代過程。

---

參考：為什么最優(yōu)解是 $(?1/\sqrt{3},?1/\sqrt{3},?1/\sqrt{3})$？

因?yàn)樵趩挝磺蚣s束下，若要最小化$x_1+x_2+x_3$，相當(dāng)于“在球面上找與(1,1,1)方向夾角最大的點(diǎn)”——也就是與 $(1,1,1)$ 反方向的單位向量。它正好是 $\frac{?1}{\sqrt{3}}(1,1,1)$。目標(biāo)值是
$$x_1+x_2+x_3=?\sqrt{3}.$$
這與幾何直覺、拉格朗日乘子法都能得到一樣的結(jié)果。

---

6. 結(jié)語

• 該代碼很好地演示了使用對數(shù)障礙(log-barrier)的內(nèi)點(diǎn)法思路：用一系列的無約束子問題（帶障礙項(xiàng)）來逼近有約束優(yōu)化，并在外層迭代中逐漸減小障礙系數(shù) (\mu)，使解貼近約束邊界的最優(yōu)解。
• 不過，在正式應(yīng)用時(shí)，往往不會(huì)用固定步長的簡單梯度下降，而會(huì)用牛頓法或線搜索，以獲得更好的數(shù)值穩(wěn)定性和收斂速度。
• 代碼本身最大的問題是 f_3d 與注釋/公式不匹配，只要改回 f_3d(x) = x(1)+x(2)+x(3) 即可與文檔保持一致，也能正確體現(xiàn)“最小化$\sum x_i$”的意圖。

在使用對數(shù)障礙（log‐barrier）形式的內(nèi)點(diǎn)法時(shí)，$\mu$（有時(shí)也記作 $t$ 或 $\nu$）通常被稱為障礙參數(shù)（barrier parameter）。它的核心作用是：

1. 控制對數(shù)障礙項(xiàng)的“強(qiáng)度”
   在障礙型目標(biāo)函數(shù)
   $$B(x,\mu )=f(x)?\mu \ln (h(x))$$
   中，(\mu) 決定了 (-\mu \ln\bigl(h(x)\bigr)) 這部分懲罰力度的大小。

   • 當(dāng) (\mu) 較大時(shí)，對(\ln(h(x)))的懲罰力度相對較小，算法對靠近約束邊界的“敏感度”不高，所以在收斂初期可以更自由地在可行域里移動(dòng)。
   • 當(dāng) (\mu) 變小時(shí)，(-\mu \ln(h(x)))會(huì)變得更尖銳，迫使解更加貼近且“貼合”可行域的邊界（若這是最優(yōu)解所在的位置）。

2. 幫助逐步逼近最優(yōu)解并保持可行
   內(nèi)點(diǎn)法的思路是：一開始用較大的 (\mu)（障礙作用較弱）來保證算法穩(wěn)步地在可行域里進(jìn)行搜索；之后在外層循環(huán)中逐步減小 (\mu)，使對數(shù)障礙項(xiàng)逐漸變得陡峭，從而把解“推”到真正需要的邊界附近并逼近最優(yōu)解。

3. 數(shù)值上的平衡

   • 如果 (\mu) 過小，一開始 (-\mu \ln(h(x))) 的勢壘就會(huì)非常強(qiáng)，導(dǎo)致很難在可行域里移動(dòng)，且容易產(chǎn)生數(shù)值不穩(wěn)定（如(\ln(h(x)))趨向負(fù)無窮）。
   • 如果 (\mu) 過大，到收斂后期也無法精確地在邊界附近找到最優(yōu)解。所以通常會(huì)有一條**“(\mu)衰減”路徑**（比如 (\mu \leftarrow \beta \mu)，(\beta<1)）來讓解逐步逼近最優(yōu)值。

簡而言之：(\mu) 是控制“障礙強(qiáng)度”的調(diào)節(jié)器。隨著 (\mu) 從大到小的不斷衰減，解會(huì)逐漸向可行域邊界靠攏并最終獲得精確的約束最優(yōu)解。

---

完整代碼

function x_history = interior_point_3d_solve()
%{使用內(nèi)點(diǎn)法(障礙函數(shù) + 固定步長梯度下降)求解:min  f(x) = x(1) + x(2) + x(3)s.t. x(1)^2 + x(2)^2 + x(3)^2 <= 1.內(nèi)點(diǎn)法: 定義障礙型目標(biāo)B(x, mu) = f(x) - mu * log( h(x) ), 其中h(x) = 1 - (x0^2 + x1^2 + x2^2).輸出:x_history: 每個(gè)迭代得到的 (x0, x1, x2).
%}% 參數(shù)mu_init = 1.0;      % 初始障礙參數(shù)mu_decay = 0.2;     % 每輪迭代后 mu 的衰減因子alpha = 0.001;       % 梯度下降步長tol = 1e-6;         % 收斂閾值max_outer = 10;     % 外層循環(huán)(更新 mu)次數(shù)max_inner = 50;     % 每次 mu 下最大迭代次數(shù)% 初始化可行解 (x0, x1, x2)，球內(nèi)，例如原點(diǎn)x = [0; 0; 0];  mu = mu_init;x_history = [];for outer = 1:max_outerfor inner = 1:max_innerg = grad_B(x, mu);  % 計(jì)算梯度x_new = x - alpha*g;% 如果越過球邊界 => h(x_new) <= 0 => x_new^2+y^2+z^2 >=1% 簡單地縮小步長再試if h_3d(x_new) <= 1e-9x_new = x - 0.1*alpha*g;end% 收斂判斷if norm(x_new - x) < tolx = x_new;break;endx = x_new;x_history = [x_history; x']; % 記錄軌跡(行向量)end% 降低 mu 使解更逼近約束邊界mu = mu * mu_decay;if mu < 1e-12break;endend% 最后再將最終點(diǎn)加入x_history = [x_history; x'];
end%% 目標(biāo)函數(shù) f(x)
function val = f_3d(x)% f(x0, x1, x2) = x0 + x1 + x2val = x(1) + x(2) + x(3);
end%% 障礙函數(shù)項(xiàng) h(x) = 1 - (x0^2 + x1^2 + x2^2)
function val = h_3d(x)val = 1.0 - (x(1)^2 + x(2)^2 + x(3)^2);
end%% 障礙型目標(biāo) B(x, mu) = f(x) - mu*ln(h(x))
function val = B_3d(x, mu)val = f_3d(x) - mu*log(h_3d(x));
end%% B(x, mu) 的梯度
function g = grad_B(x, mu)
%{B(x) = (x0 + x1 + x2) - mu * ln(1 - r^2),其中 r^2 = x0^2 + x1^2 + x2^2=> dB/dx0 = 1 + [2 mu x0 / (1 - r^2)]=> dB/dx1 = 1 + [2 mu x1 / (1 - r^2)]=> dB/dx2 = 1 + [2 mu x2 / (1 - r^2)]
%}hx = h_3d(x);r2 = x(1)^2 + x(2)^2 + x(3)^2;% hx = 1 - r2 > 0  (只要在球內(nèi))dB_dx0 = 1.0 + (2.0 * mu * x(1) / hx);dB_dx1 = 1.0 + (2.0 * mu * x(2) / hx);dB_dx2 = 1.0 + (2.0 * mu * x(3) / hx);g = [dB_dx0; dB_dx1; dB_dx2];
end%{演示如何在 3D 中用“障礙函數(shù) + 簡單梯度下降”的內(nèi)點(diǎn)法來最小化 f(x) = x0 + x1 + x2subject to x0^2 + x1^2 + x2^2 <= 1.可行域是單位球 (x0^2 + x1^2 + x2^2 <= 1)。我們會(huì)在圖中繪制球面，并用散點(diǎn)繪制迭代軌跡。
%}% 1) 調(diào)用求解函數(shù)，得到每步迭代的解 x(k)
x_history = interior_point_3d_solve();% 2) 可視化
figure('Color','w','Name','Interior-Point 3D Demo');
hold on; grid on; axis equal;  % 3D 坐標(biāo)中，最好設(shè) axis equal% 2.1 繪制單位球面（x0^2 + x1^2 + x2^2 = 1）
[Xs, Ys, Zs] = sphere(50);  
% sphere() 生成一個(gè)半徑為 1 的球面網(wǎng)格
surf(Xs, Ys, Zs, 'FaceAlpha',0.1, 'EdgeColor','none', 'FaceColor','c');
% 給球面一個(gè)半透明的青色xlabel('x_0'); ylabel('x_1'); zlabel('x_2');
title('Minimize x_0 + x_1 + x_2 subject to x_0^2 + x_1^2 + x_2^2 \le 1');% 2.2 繪制迭代軌跡
x0_hist = x_history(:,1);
x1_hist = x_history(:,2);
x2_hist = x_history(:,3);nPoints = size(x_history,1);
if nPoints > 1% 中間過程點(diǎn)用藍(lán)色散點(diǎn)plot3(x0_hist(1:end-1), x1_hist(1:end-1), x2_hist(1:end-1), ...'bo-', 'LineWidth',1.5, 'MarkerSize',4, 'MarkerFaceColor','b');
end% 最后一點(diǎn)用紅色標(biāo)記
plot3(x0_hist(end), x1_hist(end), x2_hist(end), ...'ro', 'MarkerSize',8, 'MarkerFaceColor','r');legend('Unit Sphere (Constraint)', 'Iter Process', 'Final Solution');
view(35, 25);  % 調(diào)整3D視角