如何使用ipopt进行非线性约束求目标函数最小值（NLP非线性规划）内点法（inner point method）

非线性规划，一般用matlab调用cplex和gurobi了，但这两个一般用于线性规划和二次规划

线性规划LP，二次规划（quadratic programming），如果要求更一般的非线性规划IPOT是个很好的选择，求解器很多，根据情况自己选择

非线性

具体的，这篇文章介绍的很清楚了https://blog.csdn.net/mpt0816/article/details/127638557

我这里就是再选择一个问题进行求解

ipopt的可执行程序下载下来， Releases · coin-or/Ipopt · GitHub

建立一个vs2022的工程，把include加到目录里面，把lib库都加进去，同样dll也准备好

 

就这一个主文件放入工程

编译运行即可

四个自变量，两个约束

 eval_f: 计算目标函数值，即需要最小化的目标。

 eval_grad_f: 计算目标函数的梯度。分别是4个偏导数

 eval_g: 计算约束条件的值。 n 是变量个数，m是约束条件个数，g是具体的约束函数

 eval_jac_g: 计算约束条件的雅可比矩阵（两个约束条件的一阶偏导数）

 eval_h: 计算目标函数和约束条件的二阶导数（即Hessian矩阵，二阶偏导数）。

现在使用matlab符号函数把 涉及到 用的 梯度、黑森矩阵都求一下

%

clear all
close all
clc

% 使用符号函数进行 求解梯度，黑森矩阵

syms f g1 g2
syms x1 x2 x3 x4


% 定义目标函数
f = x1 * x4 * (x1 + x2 + x3) + x3;

% 定义约束函数
g1 = x1 * x2 * x3 * x4;
g2 = x1^2 + x2^2 + x3^2 + x4^2;

% 计算目标函数的梯度和 Hessian
grad_f = gradient(f, [x1, x2, x3, x4]);
hess_f = hessian(f, [x1, x2, x3, x4]);

% 计算约束函数 g1 的梯度和 Hessian
grad_g1 = gradient(g1, [x1, x2, x3, x4]);
hess_g1 = hessian(g1, [x1, x2, x3, x4]);

% 计算约束函数 g2 的梯度和 Hessian
grad_g2 = gradient(g2, [x1, x2, x3, x4]);
hess_g2 = hessian(g2, [x1, x2, x3, x4]);

得到如下结果:

目标函数 f 的梯度：
 x1*x4 + x4*(x1 + x2 + x3)
                     x1*x4
                 x1*x4 + 1
         x1*(x1 + x2 + x3)
 
目标函数 f 的 Hessian：
[           2*x4, x4, x4, 2*x1 + x2 + x3]
[             x4,  0,  0,             x1]
[             x4,  0,  0,             x1]
[ 2*x1 + x2 + x3, x1, x1,              0]
 
约束函数 g1 的梯度：
 x2*x3*x4
 x1*x3*x4
 x1*x2*x4
 x1*x2*x3

 约束函数 g2 的梯度：
 2*x1
 2*x2
 2*x3
 2*x4 

从g1 g2看出来

   nele_jac = 8; 8个非零，两个约束条件，4个变量
   nele_hess = 10;  4*5/2=10，看其中一个hess矩阵的上三角阵

约束函数 g1 的 Hessian：
[     0, x3*x4, x2*x4, x2*x3]
[ x3*x4,     0, x1*x4, x1*x3]
[ x2*x4, x1*x4,     0, x1*x2]
[ x2*x3, x1*x3, x1*x2,     0]
 

约束函数 g2 的 Hessian：
[ 2, 0, 0, 0]
[ 0, 2, 0, 0]
[ 0, 0, 2, 0]
[ 0, 0, 0, 2]
 要替换的部分

1、eval_f 中 目标函数

2、eval_grad_f 中的梯度

   grad_f[0] = x[0] * x[3] + x[3] * (x[0] + x[1] + x[2]);
   grad_f[1] = x[0] * x[3];
   grad_f[2] = x[0] * x[3] + 1;
   grad_f[3] = x[0] * (x[0] + x[1] + x[2]);

3、eval_g 约束条件
   g[0] = x[0] * x[1] * x[2] * x[3] + my_data->g_offset[0];
   g[1] = x[0] * x[0] + x[1] * x[1] + x[2] * x[2] + x[3] * x[3] + my_data->g_offset[1];

4、eval_jac_g 约束函数的jacobi矩阵

if中 （8个）,位置是

00  01 02 03

10 11 12 13，

else 中

g1梯度，g2梯度

5、eval_h 黑森矩阵

固定抄写，4是变量个数

目标的黑森矩阵，注意走位，注意骚走位，注意下三角阵骚走位

约束的黑森

6、主函数

初始值和 上下限

约束条件的jacobi矩阵和hess矩阵的非零元素

8个=2*自变量个数

10个=自变量个数*（自变量个数+1）/2

初始值

matlab符号函求解出来的各种算式写成c有点麻烦，我这边搞了一个函数可以很方便转c

function f_str = changetoc(f)


syms x1 x2 x3 x4 %替换c语言风格
syms R %为了 R^2也能转

% f = x1^2 + x2^2 + x3^2 + x4^2;  % 示例符号函数
% f = x1^2 + x2^2 + (x1 + x2)^2 + x3^2 + x4^2;  % 示例符号函数，包含复杂表达式
% f = (r*sin(theta)*(3*cos(x1) - 4) + (x2*cos(theta)*(2*cos(x1) - 2))/n1 + (x2*sin(theta)*sin(x1))/n1)^2

% 将符号函数转换为字符串
f_str = char(f);

% 替换变量为 C 风格的数组索引 x[0], x[1], x[2], x[3]
f_str = regexprep(f_str, 'x1', 'x[0]');
f_str = regexprep(f_str, 'x2', 'x[1]');
f_str = regexprep(f_str, 'x3', 'x[2]');
f_str = regexprep(f_str, 'x4', 'x[3]');

% 定义一个集合（Cell数组）用于保存普通变量名
variables = {'x[0]','x[1]','x[2]','x[3]', 'R'};

% 
% % 示例复杂表达式
% f = (r*sin(theta)*(3*cos(x1) - 4) + (x2*cos(theta)*(2*cos(x1) - 2))/n1 + (x2*sin(theta)*sin(x1))/n1)^2 - R^2 + ...
%     (r*cos(theta) + r*sin(theta)*(6*x1 - 6*sin(x1)) + (x2*sin(theta)*(2*cos(x1) - 2))/n1 + ...
%     (x2*cos(theta)*(3*x1 - 4*sin(x1)))/n1)^2;

% 将符号函数转换为字符串
% f_str = char(f);

% 1. 替换普通变量的平方为自乘形式
for i = 1:length(variables)
    % 构建正则表达式，匹配形如 x1^2, x2^2 等
    var_pattern = strcat(variables{i}, '^2');
    % 构建替换字符串 (x1*x1), (x2*x2)
    replacement = strcat('(', variables{i}, '*', variables{i}, ')');
    % 进行替换
    f_str = strrep(f_str, var_pattern, replacement);
end

% % 找到 x[i]^2 形式的幂运算，并替换为 (x[i]*x[i])
f_str = regexprep(f_str, '(\w+\[\d+\])\^2', '$1*$1');

% 2. 替换括号表达式的平方为自乘形式
% 匹配 (xxxx)^2，替换为 (xxxx)*(xxxx)
% f_str = regexprep(f_str, '\(([^\)]+)\)\^2', '($1)*($1)');
f_str = regexprep(f_str, '\((.*?)\)\^2', '($1)*($1)');

% 输出替换后的表达式
disp(f_str);


end