Matlab 最优化求解

文章目录

@[toc]
1 线性规划
2 非线性规划
3 无约束问题matlab解法
- 3.1 符号解
- 3.2 数值解
4 求解零点与方程组的解
- 4.1 多项式
- 4.2 符号解
- 4.3 数值解
- 4.4 方程组情形
5 约束问题matlab解法
- 5.1 二次规划
- 5.2 罚函数法（外）
6 其他一些命令
- 6.1 fminbnd函数
- 6.2 fseminf函数
- 6.3 fminimax函数

matlab求解最优化函数主要包括linprog、quadprog、fminbnd、fmincon、fseminf和fminmax

1 线性规划

目标函数与约束条件都为线性函数，其Matlab标准式为
m i n f T x min \boldsymbol f^{T}\boldsymbol x minfTx

s . t . { A x ≤ b A e q ⋅ x = b e q l b ≤ b ≤ u b s.t.\left\{\begin{array}{l} \boldsymbol A\boldsymbol x \le\boldsymbol b\\ Aeq \cdot \boldsymbol x = beq\\ lb \le \boldsymbol b \le ub \end{array}\right. s.t.⎩⎨⎧Ax≤bAeq⋅x=beqlb≤b≤ub

其中 A 、 A e q \boldsymbol A、Aeq A、Aeq表示分别为不等式与等式约束系数矩阵； f , x , b , b e q , l b , u b \boldsymbol f,\boldsymbol x ,\boldsymbol b,beq,lb,ub f,x,b,beq,lb,ub为列向量。matlab求解代码为

[x,fval] = linprog(f,A,b,aeq,beq,lb,ub)

只能求出目标函数最小值，如果要求最大值，需要对目标函数转化

例1求解线性规划
max ⁡ z = 2 x 1 + 3 x 2 − 5 x 3 \max z = 2x_1+3x_2-5x_3 maxz=2x1+3x2−5x3

s . t . { x 1 + x 2 + x 3 = 7 2 x 1 − 5 x 2 + x 3 ≥ 10 x 1 + 3 x 2 + x 3 ≤ 12 x 1 , x 2 , x 3 ≥ 0 s.t. \left\{\begin{array}{l} x_1+x_2+x_3 = 7\\ 2x_1-5x_2+x_3\ge10\\ x_1+3x_2+x_3\le12\\ x_1,x_2,x_3\ge 0 \end{array}\right. s.t.⎩⎪⎪⎨⎪⎪⎧x1+x2+x3=72x1−5x2+x3≥10x1+3x2+x3≤12x1,x2,x3≥0

求解代码：

% 规划问题与matlab
%---------------------线性规划------------------
% 线性规划求解
% 例1
clc;clear;
f = [-2;-3;5]; %matlab的向量为列向量
a = [-2,5,-1; 1,3,1]; %系数矩阵
b = [-10,12]; %常数项
aeq = [1,1,1]; % 等式矩阵
beq = 7; %等式对应的常数向量
[x,y] = linprog(f,a,b,aeq,beq,zeros(3,1)) %zeros(3,1)解的下限都为正

例2 求解线性规划
max ⁡ z = 3 x 1 − x 2 − x 3 \max z = 3x_1-x_2-x_3 maxz=3x1−x2−x3

s . t . { x 1 − 2 x 2 + x 3 ≤ 11 − 4 x 1 + x 2 + 2 x 3 ≥ 3 − 2 x 1 + x 3 = 1 x 1 , x 2 , x 3 ≥ 0 s.t. \left\{\begin{array}{l} x_1-2x_2+x_3\le 11\\ -4x_1+x_2+2x_3\ge 3\\ -2x_1+x_3 = 1\\ x_1,x_2,x_3\ge 0 \end{array}\right. s.t.⎩⎪⎪⎨⎪⎪⎧x1−2x2+x3≤11−4x1+x2+2x3≥3−2x1+x3=1x1,x2,x3≥0

求解代码

clear;
f = [3;-1;-1];
a = [1,-2,1;4,-1,-2];
b = [11,-3];
aeq = [-2,0,1];
beq = 1;
[x,y] = linprog(-f,a,b,aeq,beq,zeros(3,1))

例3 （参数情形）求解线性规划
max ⁡ z = 3 x 1 − x 2 − x 3 \max z = 3x_1-x_2-x_3 maxz=3x1−x2−x3

s . t . { a x 1 − 2 x 2 + x 3 ≤ 11 − 4 x 1 + x 2 + 2 x 3 ≥ 3 − 2 x 1 + x 3 = 1 x 1 , x 2 , x 3 ≥ 0 s.t. \left\{\begin{array}{l} ax_1-2x_2+x_3\le 11\\ -4x_1+x_2+2x_3\ge 3\\ -2x_1+x_3 = 1\\ x_1,x_2,x_3\ge 0 \end{array}\right. s.t.⎩⎪⎪⎨⎪⎪⎧ax1−2x2+x3≤11−4x1+x2+2x3≥3−2x1+x3=1x1,x2,x3≥0

其中 0 < a < 10 0<a<10 0<a<10

clear;
i = 0;
hold on
while i < 10f = [3;-1;-1];a = [i*1,-2,1;4,-1,-2];b = [11,-3];aeq = [-2,0,1];beq = 1;[x,y] = linprog(-f,a,b,aeq,beq,zeros(3,1));y = -y;plot(i,y,"*r");i = i + 0.05;
end
xlabel("参数");ylabel("目标函数值")
hold off

2 非线性规划

matlab非线性规划标准式
min ⁡ f ( x ) \min f(\boldsymbol x) minf(x)

s . t { A ⋅ x ≤ b A e q ⋅ x = b e q c ( x ) ≤ 0 c e q ( x ) = 0 l b ≤ x ≤ u b s.t \left\{\begin{array}{l} \boldsymbol A\cdot \boldsymbol x \le \boldsymbol b\\ Aeq \cdot \boldsymbol x = beq\\ c(\boldsymbol x)\le 0\\ ceq(\boldsymbol x) = 0\\ lb \le \boldsymbol x \le ub \end{array}\right. s.t⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧A⋅x≤bAeq⋅x=beqc(x)≤0ceq(x)=0lb≤x≤ub

其中 c ( x ) c(\boldsymbol x) c(x)表示非线性不等式约束； c e q ( x ) ceq(\boldsymbol x) ceq(x)表示非线性等式约束。matlab一般求解代码为

[x,fval] = fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon,options)

其中nonlcon是用M文件定义的非线性向量函数 c ( x ) c(\boldsymbol x) c(x)和 c e q ( x ) ceq(\boldsymbol x) ceq(x)，options为优化参数选项。

例4 求解非线性规划

M文件1:目标函数

function f = fun1(x);
% 非线性规划求解
% 目标函数
f = sum(x.^2)+8;
end

M文件2：非线性约束条件

function [g,h] = fun2(x,a)
% 非线性规划求解
% 目标函数
g = [-x(1)^2+x(2)-x(3)^2x(1)+x(2)^2+x(3)^3 - 20]; % 非线性不等式约束
h = [-x(1)-x(2)^2 + 2x(2)+2*x(3)^2 - 3];       % 非线性等式约束
end

% 调用fmincon函数
[x,y] = fmincon("fun1",rand(3,1),[],[],[],[],zeros(3,1),[],"fun2")

3 无约束问题matlab解法

3.1 符号解

例5 求多元函数 f ( x , y ) = x 3 − y 3 + 3 x 2 + 3 y 2 − 9 x f(x,y) = x^3-y^3+3x^2+3y^2-9x f(x,y)=x3−y3+3x2+3y2−9x

% 符号解
clc,clear;
% 符号变量定义
syms x y
f = x^3 - y^3 + 3*x^2 + 3*y^2 - 9*x;
% 求梯度
df = jacobian(f);
% 求hessian矩阵
d2f = jacobian(df);
% 求驻点
[xx,yy] = solve(df);
% 转双精度类型
xx = double(xx);yy = double(yy);for i = 1:length(xx)% hessian矩阵a = subs(d2f,{x,y},{xx(i),yy(i)});% 特征值b = eig(a);f = subs(f,{x,y},{xx(i),yy(i)});if all(b > 0)fprintf("(%f,%f)是极小值点，对应极小值为% f\n",xx(i),yy(i),f);elseif all(b< 0)fprintf("(%f,%f)是极大值点，对应极大值为% f\n",xx(i),yy(i),f);elseif any(b>0) && any(b < 0)fprintf("(%f,%f)不是极值点\n",xx(i),yy(i));else fprintf("无法判断(%f,%f)是极值点\n",xx(i),yy(i));end
end

3.2 数值解

例6 求多元函数 f ( x , y ) = x 3 − y 3 + 3 x 2 + 3 y 2 − 9 x f(x,y) = x^3-y^3+3x^2+3y^2-9x f(x,y)=x3−y3+3x2+3y2−9x

% 数值解法解
clc,clear
% 定义匿名函数
f = @(x) x(1)^3 - x(2)^3 + 3*x(1)^2 + 3*x(2)^2 - 9*x(1);
g = @(x) -f(x);
% 求极小值点
[xy1,z1] = fminunc(f,rand(2,1))
% 求极大值点
[xy2,z2] = fminunc(g,rand(2,1));
xy2,z2 = -z2

利用数值方法可以加入梯度选项，编写M文件3

function [f,g] = fun3(x)
% 定义梯度
f = 100*(x(2) - x(1)^2)^2 + (1 - x(1))^2;
g = [-400*x(1)*(x(2) - x(1)^2) - 2*(1 - x(1));200*(x(2)-x(1)^2)];
end

调用函数fminunc

% 梯度法
options = optimset("GradObj","on");
[x,y] = fminunc("fun3",rand(1,2),options)

也可以加入Hessian矩阵信息，编写M文件4

function [f,g] = fun4(x)
% 定义函数
f = 100*(x(2) - x(1)^2)^2 + (1 - x(1))^2;
% 求梯度
df = [-400*x(1)*(x(2) - x(1)^2) - 2*(1 - x(1));200*(x(2)-x(1)^2)];
% 求Hessian矩阵
d2f = [-400*x(2)+1200*x(1)^2, -400*x(1)-400*x(1), 200];
end

调用fminsearch函数

% Hessian矩阵法
clc,clear
options = optimset("GradObj","on","Hessian","on");
[f,GRAD,HESS] = fminsearch("fun4",rand(1,2),options)

例7：求函数 f ( x ) = s i n ( x ) + 3 f(x) = sin(x)+3 f(x)=sin(x)+3取极小值的 x x x值

% 定义匿名函数
clc,clear
f = @(x) sin(x) + 3;
x0 = 2;
[x,y] = fminsearch(f,x0)

4 求解零点与方程组的解

4.1 多项式

例8 求解 f ( x ) = x 3 − x 2 + 2 x − 3 f(x) = x^3-x^2+2x-3 f(x)=x3−x2+2x−3的根

% 多项式系数(降幂)
xishu = [1 -1 2 -3];
x0 = roots(xishu)

4.2 符号解

% 符号解
syms x
x0 = solve(x^3 - x^2 + 2*x -3)
% 转化为小数格式,6位小数
x0 = vpa(x0,6)

4.3 数值解

% 数值解（仅求出实根）
clc,clear
y = @(x)x^3 - x^2 + 2*x -3;
x = fsolve(y,rand)

4.4 方程组情形

符号解

clc,clear
syms x y
[x,y] = solve(x^2+y-6,y^2+x-6)

数值解

f = @(x)[x(1)^2+x(2)-6,x(2)^2+x(1)-6];
% 只能求初始值需要在附近（局部）
xy = fsolve(f,rand(2,1))

5 约束问题matlab解法

5.1 二次规划

二次规划是指目标函数为二次函数但约束条件是线性的。matlab标准形式
min ⁡ 1 2 x T H x + f T x \min \frac{1}{2}\boldsymbol x^{T}\boldsymbol H\boldsymbol x+ \boldsymbol f^{T}\boldsymbol x min21xTHx+fTx

s . t . { A x ≤ b A e q ⋅ x = b e q l b ≤ b ≤ u b s.t.\left\{\begin{array}{l} \boldsymbol A\boldsymbol x \le\boldsymbol b\\ Aeq \cdot \boldsymbol x = beq\\ lb \le \boldsymbol b \le ub\end{array}\right. s.t.⎩⎨⎧Ax≤bAeq⋅x=beqlb≤b≤ub

clc,clear
h = [4,-4;-4,8];
f = [-6;-3];
a = [1,1;4,1];
b = [3;9];
[x,value] = quadprog(h,f,a,b,[],[],zeros(2,1))

5.2 罚函数法（外）

利用罚函数可以将非线性约束规划转化为求解系列无约束极值问题。也称这种方法为SUMT方法.matlab标准问题
min ⁡ f ( x ) \min f(x) minf(x)

s . t . { g i ( x ) ≤ 0 i = 1 , … , r h j ( x ) ≥ 0 j = 1 , … , s k m ( x ) = 0 m = 1 , … , t s.t.\left\{\begin{array}{l} g_i(x)\le 0 & i = 1,\dots,r\\ h_j(x)\ge 0 & j = 1,\dots ,s\\ k_m(x) =0 & m = 1,\dots ,t \end{array}\right. s.t.⎩⎨⎧gi(x)≤0hj(x)≥0km(x)=0i=1,…,rj=1,…,sm=1,…,t

设充分大 M > 0 M>0 M>0的数，构造罚函数
P ( x , M ) = f ( x ) + M ∑ i = 1 r m a x ( g i ( x ) , 0 ) − M ∑ j = 1 s m i n ( h j ( x ) , 0 ) + M ∑ m = 1 t ∣ k m ( x ) ∣ P(x,M) = f(x)+M\sum_{i=1}^rmax(g_i(x),0)-M\sum_{j=1}^smin(h_j(x),0)+M\sum_{m=1}^t|k_m(x)| P(x,M)=f(x)+Mi=1∑rmax(gi(x),0)−Mj=1∑smin(hj(x),0)+Mm=1∑t∣km(x)∣
例9求解非线性规划
min ⁡ f ( x ) = x 1 2 + x 2 2 + 8 \min f(x) = x_1^2+x_2^2+8 minf(x)=x12+x22+8

s . t { x 1 2 − x 2 ≥ 0 − x 1 − x 2 2 + 2 = 0 x 1 , x 2 ≥ 0 s.t \left\{\begin{array}{l} x_1^2-x_2\ge 0\\ -x_1-x_2^2+2 = 0\\ x_1,x_2\ge 0 \end{array}\right. s.t⎩⎨⎧x12−x2≥0−x1−x22+2=0x1,x2≥0

定义增广目标函数，编写M文件

function g = test1(x)
%罚函数构建
M = 50000;
f = x(1)^2+x(2)^2+8;
g = f-M*min(x(1),0)-M*min(x(2),0)-M*min(x(1)^2-x(2),0)+M*abs(-x(1)-x(2)^2+2);
end

调用函数fminunc

clc,clear
% 很难求出全局最优解，每次都是局部最优；每次结果存在差异
[x,y] = fminunc("test1",rand(2,1))

6 其他一些命令

6.1 fminbnd函数

用于求解单变量非线性函数在区间上的极小值，matlab格式命令

[x,fval] = fminbnd(fun,x1,x2,options)

例10 求函数 f ( x ) = ( x − 3 ) 2 − 1 , x ∈ [ 0 , 5 ] f(x) = (x-3)^2-1,x\in[0,5] f(x)=(x−3)2−1,x∈[0,5]的最小值

clc,clear
f = @(x) (x-3)^2-1;
% 在区间（0，5）上的极小值
[x,y] = fminbnd(f,0,5)

6.2 fseminf函数

用于求标准问题
min ⁡ f ( x ) s . t { A ⋅ x ≤ b A e q ⋅ x = b e q c ( x ) ≤ 0 c e q ( x ) = 0 l b ≤ x ≤ u b K i ( x , w i ) ≤ 0 , 1 ≤ i ≤ n \min f(x)\\ s.t \left\{\begin{array}{l} \boldsymbol A\cdot \boldsymbol x \le \boldsymbol b\\ Aeq \cdot \boldsymbol x = beq\\ c(\boldsymbol x)\le 0\\ ceq(\boldsymbol x) = 0\\ lb \le \boldsymbol x \le ub\\ K_i(\boldsymbol x,w_i)\le 0,1\le i\le n \end{array}\right. minf(x)s.t⎩⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎧A⋅x≤bAeq⋅x=beqc(x)≤0ceq(x)=0lb≤x≤ubKi(x,wi)≤0,1≤i≤n
例11 求函数 f ( x ) = ( x 1 − 0.5 ) 2 + ( x 2 − 0.5 ) 2 + ( x 3 − 0.5 ) 2 f(x) = (x_1-0.5)^2+(x_2-0.5)^2+(x_3-0.5)^2 f(x)=(x1−0.5)2+(x2−0.5)2+(x3−0.5)2取最小值时 x x x的值，约束条件
{ K 1 ( x , w 1 ) = s i n ( w 1 x 1 ) c o s ( w 1 x 2 ) − 1 1000 ( w 1 − 50 ) 2 − s i n ( w 1 x 3 ) − x 3 ≤ 1 K 2 ( x , w 2 ) = s i n ( w 2 x 2 ) c o s ( w 2 x 1 ) − 1 1000 ( w 2 − 50 ) 2 − s i n ( w 2 x 3 ) − x 3 ≤ 1 1 ≤ w 1 ≤ 100 , 1 ≤ w 2 ≤ 100 \left\{\begin{array}{l} K_1(x,w_1) = sin(w_1x_1)cos(w_1x_2)-\frac{1}{1000}(w_1-50)^2-sin(w_1x_3)-x_3\le 1\\ K_2(x,w_2) = sin(w_2x_2)cos(w_2x_1)-\frac{1}{1000}(w_2-50)^2-sin(w_2x_3)-x_3\le 1\\ 1\le w_1 \le 100,1\le w_2\le 100 \end{array}\right. ⎩⎨⎧K1(x,w1)=sin(w1x1)cos(w1x2)−10001(w1−50)2−sin(w1x3)−x3≤1K2(x,w2)=sin(w2x2)cos(w2x1)−10001(w2−50)2−sin(w2x3)−x3≤11≤w1≤100,1≤w2≤100
编写目标函数M文件

function f = fun5(x,s)
% 目标函数
f = sum((x - 0.5).^2);
end

编写约束函数M文件

function [c,ceq,k1,k2,s] = fun6(x,s)
% 半无穷约束
c = [];ceq = [];
if isnan(s(1,1))s = [0.2,0;0.2,0];
end
% 取样值
w1 = 1:s(1,1):100;
w2 = 1:s(2,1):100;
% 半无穷约束
k1 = sin(w1*x(1)).*cos(w1*x(2))-1/1000*(w1-50).^2-sin(w1*x(3))-x(3)-1;
k2 = sin(w2*x(2)).*cos(w2*x(1))-1/1000*(w2-50).^2-sin(w2*x(3))-x(3)-1;
% 画出版无穷约束图形
figure(2)
plot(w1,k1,"-",w2,k2,"+")
end

调用函数fseminf

clc,clear
[x,y] = fseminf(@fun5,rand(3,1),2,@fun6)

6.3 fminimax函数

用于求解标准问题
min ⁡ x max ⁡ i F i ( x ) s . t { A ⋅ x ≤ b A e q ⋅ x = b e q c ( x ) ≤ 0 c e q ( x ) = 0 l b ≤ x ≤ u b \min_x\max_i F_i(x)\\ s.t \left\{\begin{array}{l} \boldsymbol A\cdot \boldsymbol x \le \boldsymbol b\\ Aeq \cdot \boldsymbol x = beq\\ c(\boldsymbol x)\le 0\\ ceq(\boldsymbol x) = 0\\ lb \le \boldsymbol x \le ub \end{array}\right. xminimaxFi(x)s.t⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧A⋅x≤bAeq⋅x=beqc(x)≤0ceq(x)=0lb≤x≤ub
matlab命令为

[x,fval] = fminimax(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon,options)

例11求解函数族 { f 1 ( x ) , f 2 ( x ) , f 3 ( x ) , f 4 ( x ) , f 5 ( x ) } \{f_1(x),f_2(x),f_3(x),f_4(x),f_5(x)\} {f1(x),f2(x),f3(x),f4(x),f5(x)}取极小-极大值时 X X X的值，其中
{ f 1 ( x ) = 2 x 1 2 + x 2 2 − 48 x 1 − 40 x 2 + 304 f 2 ( x ) = − x 1 2 − 3 x 2 2 f 3 ( x ) = x 1 + 3 x 2 − 18 f 4 ( x ) = − x 1 − x 2 f 5 ( x ) = x 1 + x 2 − 8 \left\{\begin{array}{l} f_1(x) = 2x_1^2+x_2^2-48x_1-40x_2+304\\ f_2(x) = -x_1^2-3x_2^2\\ f_3(x) = x_1+3x_2-18\\ f_4(x) = -x_1-x_2\\ f_5(x) = x_1+x_2-8 \end{array}\right. ⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧f1(x)=2x12+x22−48x1−40x2+304f2(x)=−x12−3x22f3(x)=x1+3x2−18f4(x)=−x1−x2f5(x)=x1+x2−8
定义向量函数M文件

function f = fun7(x)
% 目标函数
f = [2*x(1)^2+x(2)^2-48*x(1)-40*x(2)+304-x(1)^2 - 3*x(2)^2x(1)+3*x(2)-18-x(1)-x(2)x(1)+x(2)-8];
end

调用函数fminimax

clc,clear
[x,y] = fminimax(@fun7,rand(2,1))

-END-

参考文献

司守奎等.数学建模算法与应用[M].国防工业出版社