咨询求二维函数f(a,b)最大值的高效数值方法及C实现

你的函数计算耗时极长且属于黑箱场景（未提及可获取导数信息），优先推荐Nelder-Mead单纯形法——它无需计算函数梯度/导数，仅通过函数值比较迭代优化，能有效减少函数调用次数，非常适配这类需求。另外也可采用网格搜索+一维优化的嵌套策略，但Nelder-Mead的收敛效率通常更高。

核心方法说明 #

Nelder-Mead单纯形法 #

该方法在二维空间中维护由3个顶点组成的单纯形，通过比较顶点函数值，执行反射、扩展、收缩等操作调整单纯形的形状与位置，逐步逼近极值点。每轮迭代的函数调用次数有限，无需导数支撑，是黑箱函数优化的首选方案之一。

嵌套一维优化 #

若能接受稍多调用量，可采用交替迭代的嵌套逻辑：

• 固定参数a，用黄金分割法（比牛顿法更适合无导数场景）找到使f(a,b)最大的b
• 固定找到的b，再用一维优化寻找最优a
• 重复上述步骤直至收敛
  该方法逻辑简单，但收敛速度通常不如Nelder-Mead。

C语言实现示例（Nelder-Mead求最大值） #

注：通过将目标函数取反，把求最大值转化为Nelder-Mead擅长的最小值求解问题。

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

// 替换为你的目标函数f(a,b)，示例函数为f(a,b) = -(a² + b²)，最大值在(0,0)处
double target_function(double a, double b) {
    return -(a*a + b*b); // 取反后求最小值等价于原函数求最大值
}

// 计算单纯形所有顶点的函数值
void evaluate_simplex(double simplex[3][2], double values[3]) {
    for (int i = 0; i < 3; i++) {
        values[i] = target_function(simplex[i][0], simplex[i][1]);
    }
}

// 排序找到单纯形中函数值最小、次小、最大的顶点索引
void sort_simplex(double values[3], int *idx_low, int *idx_mid, int *idx_high) {
    *idx_low = 0;
    *idx_mid = 1;
    *idx_high = 2;

    if (values[1] < values[0]) {
        *idx_low = 1;
        *idx_mid = 0;
    }
    if (values[2] < values[*idx_low]) {
        *idx_high = *idx_mid;
        *idx_mid = *idx_low;
        *idx_low = 2;
    } else if (values[2] < values[*idx_mid]) {
        *idx_high = *idx_mid;
        *idx_mid = 2;
    }
}

// Nelder-Mead优化主逻辑
void nelder_mead(double start_a, double start_b, double tol, int max_iter, double *opt_a, double *opt_b) {
    // 初始化单纯形（可根据参数范围调整初始步长）
    double simplex[3][2] = {
        {start_a, start_b},
        {start_a + 0.1, start_b},
        {start_a, start_b + 0.1}
    };
    double values[3];
    evaluate_simplex(simplex, values);

    int iter = 0;
    while (iter < max_iter) {
        int idx_low, idx_mid, idx_high;
        sort_simplex(values, &idx_low, &idx_mid, &idx_high);

        // 收敛判断：顶点函数值极差小于阈值
        double range = values[idx_high] - values[idx_low];
        if (range < tol) break;

        // 计算单纯形中心（排除最大值顶点）
        double centroid[2] = {0};
        for (int i = 0; i < 3; i++) {
            if (i != idx_high) {
                centroid[0] += simplex[i][0];
                centroid[1] += simplex[i][1];
            }
        }
        centroid[0] /= 2;
        centroid[1] /= 2;

        // 反射操作
        double reflect[2] = {2*centroid[0]-simplex[idx_high][0], 2*centroid[1]-simplex[idx_high][1]};
        double val_reflect = target_function(reflect[0], reflect[1]);

        if (val_reflect < values[idx_low]) {
            // 扩展操作
            double expand[2] = {centroid[0]+2*(reflect[0]-centroid[0]), centroid[1]+2*(reflect[1]-centroid[1])};
            double val_expand = target_function(expand[0], expand[1]);
            if (val_expand < val_reflect) {
                simplex[idx_high][0] = expand[0];
                simplex[idx_high][1] = expand[1];
                values[idx_high] = val_expand;
            } else {
                simplex[idx_high][0] = reflect[0];
                simplex[idx_high][1] = reflect[1];
                values[idx_high] = val_reflect;
            }
        } else if (val_reflect < values[idx_mid]) {
            simplex[idx_high][0] = reflect[0];
            simplex[idx_high][1] = reflect[1];
            values[idx_high] = val_reflect;
        } else {
            if (val_reflect < values[idx_high]) {
                // 外部收缩
                double contract_out[2] = {centroid[0]+0.5*(reflect[0]-centroid[0]), centroid[1]+0.5*(reflect[1]-centroid[1])};
                double val_contract_out = target_function(contract_out[0], contract_out[1]);
                if (val_contract_out < val_reflect) {
                    simplex[idx_high][0] = contract_out[0];
                    simplex[idx_high][1] = contract_out[1];
                    values[idx_high] = val_contract_out;
                } else {
                    // 内部收缩
                    for (int i = 0; i < 3; i++) {
                        if (i != idx_low) {
                            simplex[i][0] = simplex[idx_low][0] + 0.5*(simplex[i][0]-simplex[idx_low][0]);
                            simplex[i][1] = simplex[idx_low][1] + 0.5*(simplex[i][1]-simplex[idx_low][1]);
                        }
                    }
                    evaluate_simplex(simplex, values);
                }
            } else {
                // 内部收缩
                for (int i = 0; i < 3; i++) {
                    if (i != idx_low) {
                        simplex[i][0] = simplex[idx_low][0] + 0.5*(simplex[i][0]-simplex[idx_low][0]);
                        simplex[i][1] = simplex[idx_low][1] + 0.5*(simplex[i][1]-simplex[idx_low][1]);
                    }
                }
                evaluate_simplex(simplex, values);
            }
        }
        iter++;
    }

    // 获取最优顶点（对应原函数最大值）
    int idx_low;
    sort_simplex(values, &idx_low, NULL, NULL);
    *opt_a = simplex[idx_low][0];
    *opt_b = simplex[idx_low][1];
}

int main() {
    double start_a = 1.0;
    double start_b = 1.0;
    double tolerance = 1e-6;
    int max_iterations = 100;
    double opt_a, opt_b;

    nelder_mead(start_a, start_b, tolerance, max_iterations, &opt_a, &opt_b);

    printf("最优参数：a = %.6f, b = %.6f\n", opt_a, opt_b);
    printf("原函数最大值：%.6f\n", -target_function(opt_a, opt_b));

    return 0;
}

注意事项 #

1. 若能获取目标函数的偏导数，可考虑共轭梯度法或二维牛顿法，但导数计算若同样耗时，效率反而不如Nelder-Mead。
2. 初始单纯形的步长需根据你的参数范围调整，示例中0.1的步长可按需放大或缩小。
3. 收敛阈值tol和最大迭代次数max_iter可根据精度需求灵活调整。

内容的提问来源于stack exchange，提问作者dargaud