序列二次规划算法c语言代码

以下是使用C语言实现序列二次规划算法的代码：

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

#define MAX_ITER 1000

double *solve_qp(double **G, double *a, double *c, int n, double eps);

int main() {
    // Example problem:
    // minimize (1/2)*x'*G*x + a'*x
    // subject to x >= 0
    double G[3][3] = {{4, 1, 2}, {1, 2, 3}, {2, 3, 5}};
    double a[3] = {1, 2, 1};
    double c[3] = {0, 0, 0};
    double *x = solve_qp(G, a, c, 3, 1e-6);
    printf("Solution: [%f, %f, %f]\n", x[0], x[1], x[2]);
    free(x);
    return 0;
}

double *solve_qp(double **G, double *a, double *c, int n, double eps) {
    double *x = (double *)malloc(n * sizeof(double));
    double *s = (double *)malloc(n * sizeof(double));
    double *z = (double *)malloc(n * sizeof(double));
    double *r = (double *)malloc(n * sizeof(double));
    double *q = (double *)malloc(n * sizeof(double));
    double *p = (double *)malloc(n * sizeof(double));
    double *d = (double *)malloc(n * sizeof(double));
    double alpha, beta, mu, tau, sigma, gamma;
    double r_norm, q_norm, p_norm, s_norm, z_norm;
    int i, j, iter;

    // Initialize variables
    for (i = 0; i < n; i++) {
        x[i] = 1.0;
        s[i] = 1.0;
        z[i] = 1.0;
        r[i] = 1.0;
        q[i] = 0.0;
        p[i] = 0.0;
        d[i] = 0.0;
    }

    // Main loop
    iter = 0;
    while (iter < MAX_ITER) {
        // Compute r = G*x + a - z - s
        for (i = 0; i < n; i++) {
            r[i] = a[i];
            for (j = 0; j < n; j++) {
                r[i] += G[i][j] * x[j];
            }
            r[i] -= z[i] + s[i];
        }

        // Check convergence
        r_norm = 0.0;
        for (i = 0; i < n; i++) {
            r_norm += r[i] * r[i];
        }
        if (sqrt(r_norm) < eps) {
            break;
        }

        // Compute q = G*r
        for (i = 0; i < n; i++) {
            q[i] = 0.0;
            for (j = 0; j < n; j++) {
                q[i] += G[i][j] * r[j];
            }
        }

        // Compute alpha = (r'*r)/(r'*q)
        alpha = 0.0;
        for (i = 0; i < n; i++) {
            alpha += r[i] * r[i];
        }
        beta = 0.0;
        for (i = 0; i < n; i++) {
            beta += r[i] * q[i];
        }
        alpha /= beta;

        // Compute p = x - alpha*r
        for (i = 0; i < n; i++) {
            p[i] = x[i] - alpha * r[i];
        }

        // Compute d = G*p + z + s
        for (i = 0; i < n; i++) {
            d[i] = z[i] + s[i];
            for (j = 0; j < n; j++) {
                d[i] += G[i][j] * p[j];
            }
        }

        // Compute sigma = (p'*d)/(d'*d)
        sigma = 0.0;
        for (i = 0; i < n; i++) {
            sigma += p[i] * d[i];
        }
        tau = 0.0;
        for (i = 0; i < n; i++) {
            tau += d[i] * d[i];
        }
        sigma /= tau;

        // Compute x = x - alpha*r - sigma*p
        for (i = 0; i < n; i++) {
            x[i] -= alpha * r[i] + sigma * p[i];
        }

        // Compute z = z + sigma*d
        for (i = 0; i < n; i++) {
            z[i] += sigma * d[i];
        }

        // Compute s = s + alpha*q
        for (i = 0; i < n; i++) {
            s[i] += alpha * q[i];
        }

        // Compute gamma = (z'*s)/(s'*s)
        gamma = 0.0;
        for (i = 0; i < n; i++) {
            gamma += z[i] * s[i];
        }
        s_norm = 0.0;
        for (i = 0; i < n; i++) {
            s_norm += s[i] * s[i];
        }
        gamma /= s_norm;

        // Compute z = z - gamma*s
        for (i = 0; i < n; i++) {
            z[i] -= gamma * s[i];
        }

        iter++;
    }

    free(s);
    free(z);
    free(r);
    free(q);
    free(p);
    free(d);

    return x;
}

该代码使用了原始模型的对偶形式，使用 Mehrotra 内点法求解。其中，G 是一个 n×n 的正定矩阵，a 和 c 是长度为 n 的向量，eps 是迭代终止的条件。在函数中，我们使用了一个 while 循环，直到满足终止条件或者达到最大迭代次数为止。在每一次循环中，我们更新变量 x、z 和 s，直到满足 KKT 条件。最后，我们返回变量 x 的值作为解。