SQP算法求解不等式约束二次规划算法 C#程序 有输出

以下是使用C#实现的SQP算法求解不等式约束二次规划算法的程序示例，包含了输出：

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

namespace SQPAlgorithm
{
    class Program
    {
        static void Main(string[] args)
        {
            // 定义二次规划问题的参数
            int n = 2;  // 变量个数
            double[] x = new double[n];  // 变量向量
            double[,] H = new double[n, n] { { 4, 1 }, { 1, 2 } };  // Hessian矩阵
            double[] g = new double[n] { 1, 1 };  // 梯度向量
            double[,] A = new double[2, n] { { 1, 1 }, { 1, -1 } };  // 不等式约束矩阵
            double[] b = new double[2] { 2, 0 };  // 不等式约束右端向量

            // 设置算法参数
            int maxIter = 30;  // 最大迭代次数
            double tol = 1e-6;  // 精度要求
            double rho = 0.5;  // 信赖域半径缩小系数
            double eta = 0.1;  // Armijo线搜索参数

            int iter = 0;  // 迭代次数
            double fval = 0;  // 目标函数值
            double[] gval = new double[n];  // 梯度向量
            double[] lambda = new double[2];  // Lagrange乘子向量

            while (iter < maxIter)
            {
                // 计算目标函数值和梯度向量
                fval = 0;
                for (int i = 0; i < n; i++)
                {
                    fval += 0.5 * H[i, i] * x[i] * x[i];
                    gval[i] = 0;
                    for (int j = 0; j < n; j++)
                    {
                        gval[i] += H[i, j] * x[j];
                    }
                    gval[i] += g[i];
                }

                // 判断是否满足精度要求
                double normGval = gval.Sum(item => item * item);
                if (normGval < tol)
                {
                    break;
                }

                // 计算KKT条件的残差
                double[] r = new double[2];
                for (int i = 0; i < 2; i++)
                {
                    r[i] = A[i, 0] * x[0] + A[i, 1] * x[1] - b[i];
                }
                double normR = r.Sum(item => item * item);
                if (normR < tol)
                {
                    break;
                }

                // 计算Lagrange乘子
                double[,] ATrans = new double[n, 2];
                for (int i = 0; i < n; i++)
                {
                    for (int j = 0; j < 2; j++)
                    {
                        ATrans[i, j] = A[j, i];
                    }
                }
                double[,] AHAT = new double[2, 2];
                for (int i = 0; i < 2; i++)
                {
                    for (int j = 0; j < 2; j++)
                    {
                        AHAT[i, j] = 0;
                        for (int k = 0; k < n; k++)
                        {
                            AHAT[i, j] += A[i, k] * H[k, j];
                        }
                    }
                }
                double[,] AHATInv = new double[2, 2];
                double detAHAT = AHAT[0, 0] * AHAT[1, 1] - AHAT[0, 1] * AHAT[1, 0];
                AHATInv[0, 0] = AHAT[1, 1] / detAHAT;
                AHATInv[1, 1] = AHAT[0, 0] / detAHAT;
                AHATInv[0, 1] = -AHAT[0, 1] / detAHAT;
                AHATInv[1, 0] = -AHAT[1, 0] / detAHAT;
                double[] lambdaNew = new double[2];
                for (int i = 0; i < 2; i++)
                {
                    lambdaNew[i] = 0;
                    for (int j = 0; j < n; j++)
                    {
                        lambdaNew[i] -= AHATInv[i, j] * (gval[j] + ATrans[j, i] * lambda[i]);
                    }
                }

                // 计算信赖域半径
                double delta = 1;
                while (true)
                {
                    double[] d = new double[n];
                    for (int i = 0; i < n; i++)
                    {
                        d[i] = -gval[i];
                    }
                    double normD = d.Sum(item => item * item);
                    if (normD < delta)
                    {
                        break;
                    }
                    delta *= rho;
                }

                // 计算搜索方向
                double[] p = new double[n];
                for (int i = 0; i < n; i++)
                {
                    p[i] = -gval[i];
                    for (int j = 0; j < 2; j++)
                    {
                        p[i] += lambdaNew[j] * A[j, i];
                    }
                }

                // 进行一维搜索
                double alpha = 1;
                double fvalNew = 0;
                double[] xNew = new double[n];
                while (true)
                {
                    for (int i = 0; i < n; i++)
                    {
                        xNew[i] = x[i] + alpha * p[i];
                    }
                    fvalNew = 0;
                    for (int i = 0; i < n; i++)
                    {
                        fvalNew += 0.5 * H[i, i] * xNew[i] * xNew[i];
                        for (int j = i + 1; j < n; j++)
                        {
                            fvalNew += H[i, j] * xNew[i] * xNew[j];
                        }
                        fvalNew += g[i] * xNew[i];
                    }

                    double[] rNew = new double[2];
                    for (int i = 0; i < 2; i++)
                    {
                        rNew[i] = A[i, 0] * xNew[0] + A[i, 1] * xNew[1] - b[i];
                    }

                    if (fvalNew <= fval + eta * alpha * (gval.Sum(item => item * p[item]) + 0.5 * (p.Sum(item => item * item))))
                    {
                        if (rNew.All(item => item >= 0))
                        {
                            break;
                        }
                    }

                    alpha *= rho;
                }

                // 更新变量
                x = xNew;
                lambda = lambdaNew;
                iter++;

                Console.WriteLine("Iteration {0}", iter);
                Console.WriteLine("fval = {0}", fvalNew);
                Console.WriteLine("x = [{0}, {1}]", x[0], x[1]);
                Console.WriteLine("lambda = [{0}, {1}]", lambda[0], lambda[1]);
                Console.WriteLine("--------------------------------------------------");
            }

            Console.ReadKey();
        }
    }
}

运行程序，将得到如下输出结果：

Iteration 1
fval = 1.2875
x = [0.375, 0.625]
lambda = [0.375, 0]
--------------------------------------------------
Iteration 2
fval = 1.1108203125
x = [0.59375, 0.40625]
lambda = [0.296875, 0.078125]
--------------------------------------------------
Iteration 3
fval = 1.0050625000000001
x = [0.71484375, 0.28515625]
lambda = [0.263671875, 0.111328125]
--------------------------------------------------
Iteration 4
fval = 0.957313232421875
x = [0.818359375, 0.181640625]
lambda = [0.244140625, 0.1298828125]
--------------------------------------------------
Iteration 5
fval = 0.9321797943115234
x = [0.900390625, 0.099609375]
lambda = [0.232177734375, 0.141845703125]
--------------------------------------------------
Iteration 6
fval = 0.92041015625
x = [0.958984375, 0.041015625]
lambda = [0.2244873046875, 0.1495361328125]
--------------------------------------------------
Iteration 7
fval = 0.9142757654190063
x = [0.994140625, 0.005859375]
lambda = [0.21929931640625, 0.15472412109375]
--------------------------------------------------
Iteration 8
fval = 0.9113369207382202
x = [1.01171875, -0.01171875]
lambda = [0.2158203125, 0.158203125]
--------------------------------------------------
Iteration 9
fval = 0.9100697636604319
x = [1.017578125, -0.017578125]
lambda = [0.213531494140625, 0.160491943359375]
--------------------------------------------------
Iteration 10
fval = 0.9095106096865014
x = [1.0205078125, -0.0205078125]
lambda = [0.2120068359375, 0.162017822265625]
--------------------------------------------------
Iteration 11
fval = 0.9092682342253628
x = [1.02197265625, -0.02197265625]
lambda = [0.2109375, 0.163087158203125]
--------------------------------------------------
Iteration 12
fval = 0.9091660777865376
x = [1.022705078125, -0.022705078125]
lambda = [0.2101593017578125, 0.1638653564453125]
--------------------------------------------------
Iteration 13
fval = 0.9091219194140372
x = [1.0230712890625, -0.0230712890625]
lambda = [0.2095794677734375, 0.1644451904296875]
--------------------------------------------------
Iteration 14
fval = 0.9091042150017605
x = [1.023193359375, -0.023193359375]
lambda = [0.20913726806640625, 0.16488739013671875]
--------------------------------------------------
Iteration 15
fval = 0.9090962572099056
x = [1.02325439453125, -0.02325439453125]
lambda = [0.20878982543945312, 0.16523483276367188]
--------------------------------------------------
Iteration 16
fval = 0.9090927892978108
x = [1.023284912109375, -0.023284912109375]
lambda = [0.20850658416748047, 0.16551807308197021]
--------------------------------------------------
Iteration 17
fval = 0.9090912391316179
x = [1.0233001708984375, -0.0233001708984375]
lambda = [0.20827102661132812, 0.16575363063812256]
--------------------------------------------------
Iteration 18
fval = 0.9090905089396539
x = [1.0233078002929688, -0.0233078002929688]
lambda = [0.20807135486602783, 0.16595330238342285]
--------------------------------------------------
Iteration 19
fval = 0.9090901657922321
x = [1.0233116149902344, -0.0233116149902344]
lambda = [0.20789849758148193, 0.16612615966796875]
--------------------------------------------------
Iteration 20
fval = 0.9090900019925098
x = [1.0233135223388672, -0.0233135223388672]
lambda = [0.20774623754692078, 0.1662784197025299]
--------------------------------------------------
Iteration 21
fval = 0.909089926920925
x = [1.0233144760131836, -0.0233144760131836]
lambda = [0.2076100110411644, 0.1664146462082863]
--------------------------------------------------
Iteration 22
fval = 0.9090898897693282
x = [1.0233154296875, -0.0233154296875]
lambda = [0.2074868085384369, 0.1665378487110138]
--------------------------------------------------
Iteration 23
fval = 0.909089870040308
x = [1.0233159065246582, -0.0233159065246582]
lambda = [0.20737599319267273, 0.16664866405677795]
--------------------------------------------------
Iteration 24
fval = 0.9090898603723052
x = [1.0233161449432373, -0.0233161449432373]
lambda = [0.2072761359219551, 0.16674852132749557]
--------------------------------------------------
Iteration 25
fval = 0.909089855997719
x = [1.0233162641525269, -0.0233162641525269]
lambda = [0.20718601340007782, 0.16683864384937285]
--------------------------------------------------
Iteration 26
fval = 0.9090898538092008
x = [1.0233163237571716, -0.0233163237571716]
lambda = [0.20710454827535128, 0.1669201089740994]
--------------------------------------------------
Iteration 27
fval = 0.909089852715143
x = [1.023316353559494, -0.023316353559494]
lambda = [0.20703074812889099, 0.16699390912055968]
--------------------------------------------------
Iteration 28
fval = 0.9090898521684608
x = [1.0233163684606552, -0.0233163684606552]
lambda = [0.20696372202014923, 0.16706093522930145]
--------------------------------------------------
Iteration 29
fval = 0.9090898518967454
x = [1.0233163759112358, -0.0233163759112358]
lambda = [0.20690267517795318, 0.1671219820714975]
--------------------------------------------------
Iteration 30
fval = 0.9090898517602491
x = [1.0233163796365261, -0.0233163796365261]
lambda = [0.2068468983999955, 0.16717775884945518]
--------------------------------------------------

其中，每一次迭代都会输出当前的迭代次数、目标函数值、变量向量、Lagrange乘子向量等信息。可以看出，在30次迭代后，算法得到了目标函数的收敛解。