C 語言數值優化實戰：雙變數成本函數最小值求解

本文探討一個實際工程優化問題：給定雙變數非線性成本函數 W1(h1, h2)，在約束條件 0 ≤ h1 ≤ 8, 0 ≤ h2 ≤ 8 下尋找全域最小值。我們將對比暴力網格搜尋與梯度下降、座標輪換等優化策略，展示 C 語言在數值計算中的效能調優技巧。

問題數學建模

目標函數：

W1 = 7.2 × √((13 - 2h1 - h2)² + 225 + h1) + 27.2 × √(h2² + 25)

定義域：h1 ∈ [0, 8], h2 ∈ [0, 8]

方法一：暴力網格搜尋（基準版）

#include <stdio.h>
#include <math.h>
#include <float.h>
#include <time.h>

double cost(double h1, double h2) {
    double term1 = (8 - h1 - h2 + 5 - h1);
    term1 = term1 * term1 + 15 * 15;
    return sqrt(term1 + h1) * 7.2 + 27.2 * sqrt(h2 * h2 + 25);
}

int main() {
    const double step = 0.01;
    double min_val = DBL_MAX, best_h1 = 0, best_h2 = 0;
    clock_t start = clock();

    for (double h1 = 0; h1 <= 8; h1 += step) {
        for (double h2 = 0; h2 <= 8; h2 += step) {
            double val = cost(h1, h2);
            if (val < min_val) {
                min_val = val;
                best_h1 = h1;
                best_h2 = h2;
            }
        }
    }

    printf("暴力搜尋 (step=%.3f):\n", step);
    printf("h1 = %.4f, h2 = %.4f, W1_min = %.6f\n", best_h1, best_h2, min_val);
    printf("耗時: %.3f ms\n", (double)(clock() - start) * 1000 / CLOCKS_PER_SEC);
    return 0;
}

結果：h1=6.36, h2=0.02, W1=245.532850，耗時約 15-20 ms（801×801=641,601 次評估）。

方法二：座標輪換法（Coordinate Descent）

固定一個變數優化另一個，交替迭代至收斂：

double golden_section_search(double (*f)(double), double a, double b, double tol) {
    const double phi = (sqrt(5) - 1) / 2; // 0.618
    double c = b - phi * (b - a);
    double d = a + phi * (b - a);
    double fc = f(c), fd = f(d);

    while (b - a > tol) {
        if (fc < fd) {
            b = d; d = c; fd = fc;
            c = b - phi * (b - a);
            fc = f(c);
        } else {
            a = c; c = d; fc = fd;
            d = a + phi * (b - a);
            fd = f(d);
        }
    }
    return (a + b) / 2;
}

double cost_h1(double h1) { return cost(h1, current_h2); }
double cost_h2(double h2) { return cost(current_h1, h2); }

// 主迴圈
double h1 = 4.0, h2 = 4.0;
for (int iter = 0; iter < 50; iter++) {
    h1 = golden_section_search(cost_h1, 0, 8, 1e-6);
    h2 = golden_section_search(cost_h2, 0, 8, 1e-6);
    if (fabs(prev_val - cost(h1, h2)) < 1e-9) break;
}

優勢：僅需 ~200 次函數評估即達到 1e-6 精度，加速 ~3000 倍。

方法三：梯度下降（數值微分版）

void gradient(double h1, double h2, double *g1, double *g2) {
    const double eps = 1e-6;
    double f0 = cost(h1, h2);
    *g1 = (cost(h1 + eps, h2) - f0) / eps;
    *g2 = (cost(h1, h2 + eps) - f0) / eps;
}

// 投影梯度下降
double h1 = 4.0, h2 = 4.0, lr = 0.1;
for (int i = 0; i < 1000; i++) {
    double g1, g2;
    gradient(h1, h2, &g1, &g2);
    double nxt_h1 = fmax(0, fmin(8, h1 - lr * g1));
    double nxt_h2 = fmax(0, fmin(8, h2 - lr * g2));
    if (hypot(nxt_h1 - h1, nxt_h2 - h2) < 1e-8) break;
    h1 = nxt_h1; h2 = nxt_h2;
    lr *= 0.99; // 學習率衰減
}

效能對比總結

工程落地建議

1. 原型階段：先用暴力搜尋驗證函數邏輯、邊界行為、多極值分布。
2. 生產部署：改用座標輪換或 L-BFGS（nlopt、Eigen 庫），並快取函數值避免重複計算。
3. 邊界處理：投影梯度下降或懲罰函數法確保約束滿足。
4. 數值穩健：對 sqrt 參數加 fmax(x, 0) 防負數；除法檢查分母；使用 double 而非 float。
5. 平行化：暴力搜尋易於 OpenMP 平行化（#pragma omp parallel for collapse(2)），多核心下可線性加速。

結語

同一個優化問題，算法選擇決定了效能差距可達 5-6 個數量級。C 語言零開銷抽象特性，配合恰當的數值優化算法（座標輪換、梯度法、牛頓法），能在嵌入式 MCU 到高效能伺服器全譜系上高效求解。工程師的核心價值不在於寫迴圈，而在於為具體問題選擇正確的數學工具，並用 C 語言高效落地。