算法1—线性规划和非线性规划

线性规划

1.标准线性规划

• Matlab 中规定线性规划的标准形式为

其中 c 、x、b、beq、lb、ub都为列向量， A 、 Aeq 为矩阵。
Matlab中求解线性规划的命令为：

[x,fval] = linprog(c,A,b,Aeq,beq,lb,ub)
c为价值向量， x为返回变量的取值，fval为返回目标函数的最优值，A和b对应线性不等式约束，Aeq和beq对应是线性等式约束，lb和ub是变量的下界和上界向量。
若模型中没有线性不等式或者是线性等式，则在linprog()中对应参数写为空的列表（[]）

• 将线性规划转化为Matlab标准型
  线性规划为：’ 表示转置T
  max c’x
  Ax ≥ b
  Matlab标准型为：
  min -c’x
  -Ax ≤ -b
• 例题

例题1
man Z = 3X1 -X2 - X3
s.t：
X1 - 2X2 + X3 ≤ 11
-4X1 + X2 + 2X3 ≥3
-2X1 + X3 = 1
X1,X2.X3 ≥ 0

代码为：
c = [3;-1;-1]; % 列变量系数
a = [1,-2,1;4,-1,-2];
b = [11;-3];
Aeq = [-2,0,1];
beq = [1];
[x,y] = linprog(-c,a,b,Aeq,beq,zeros(3,1))
y=-y % 换算到目标函数极大化（原是求max，在标准型下求min,在最终结果基础上去相反数）

求得：X1 = 4,X2 = 1,X3 = 9, Z = 2

• zeros(M,N)就是生成M行N列的零矩阵
• ones(M,N)就是生成M行N列的元素为1的矩阵
• 单位矩阵（eye(m,n)）就是生成斜对角都为1的矩阵

2.非线性规划问题转化为线性规划

• 例题

做变量变换ui=(Xi + |Xi|)/2，vi=(|Xi| - Xi)/2，i取1,2,3,4再把新变量重新排序成一维向量y=[u,v]'=[u1;u2;u3;u4;v1;v2;v3;v4]，从而变换为线性规划模型。

例题2
min Z = |1|+2|X2|+3|X3|+4|X4|
s.t：
X1 - X2 - X3 + X4 ≤ -2
X1 - X2 + X3 -3X4 ≥ -1
X1 - X2 - 2X3 + 3X4 ≤ -1/2

代码
clc,clear
c = [1,2,3,4] % 可以写成c = 1:4
c = [c;c] % 构造价值列向量
a = [1,-1,-1,1;1,-1,1,-3;1,-1,-2,3]
a = [a, -a] % 构造变换后新的系数矩阵
b = [-2;-1;-1/2]
[y,Z] = linprog(c,a,b,[],[],zeros(8,1)) % 问题中没有等式约束，对应的矩阵为空矩阵
x = y(1:4) - y(5:end) % 变换到原问题的解，Z=u-v

输出为：X1 = -2，X2=X3=X4=0，最优值Z=2

非线性规划

1.非线性规划的matlab解法

Matlab中非线性规划的数学模型为：

c(x),ceq(x)为非线性向量函数。
Matlab命令为：

[x, fval] = fmin(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon,options)
fun是用M文件定义的函数f(x)；nonlcon是用M文件定义的非线性向量函数c(x)、ceq(x)；option定义优化参数可以默认Matlab设置。

• 例题1

% 定义目标函数fun1
function f = fun1(x)
f=sum(x.^2)+8;
% 定义非线性约束条件函数fun2
function[g,h]=fun2(x)
g=[-x(1)^2+x(2)-x(3)^2
    x(1)^2+x(2)^2+x(3)^3-20];% 非线性不等式约束
h=[-x(1)-x(2)^2+2
    x(2)+2*x(3)^2-3];% 非线性等式约束
% 主程序文件example1
[x,y]=fmincon('fun1',rand(3,1),[],[],[],[],zeros(3,1),[],'fun2')
%（所求目标函数，x的初值，没有线性约束所以A=[],b=[],Aeq=[],beq=[]，x的取值下限三个0，无上限，约束条件）

输出为：

x1=0.5522
x2=1.2033
x3=0.9478
最小值y=10.6511

2.无约束极值的Matlab解法

2.1无约束极值问题的数值解

在 Matlab 工具箱中，用于求解无约束极值问题的函数有 fminunc 和 fminsearch。
求函数的极小值：

Matlab 中 fminunc 的基本命令为：

[x, fval] = fminunc(fun,x0,options)
返回的x为极小值点，fval返回的是极小值。fun 是一个 M 文件，当 fun只有一个返回值时，它的返回值是函数 f (x)； 当 fun有两个返回值时，它的第二个返回值是 f (x)的梯度向量；当 fun 有三个返回值时，它的第三个返回值是 f (x)的二阶导数阵（Hessian 阵）。

求多元函数的极小值用fminsearch: [x, fval] = fminsearch(fun,x0,options)

• 例题2

% 求函数 f (x) = 100(x2 − x12 )^2 + (1− x1 )^2 的最小值
编写目标函数fun3
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)]; % 返回的是梯度向量
编写主程序文件example21
options = optimset('GradObj','on');
[x,y]=fminunc('fun3',rand(1,2),options)

输出为：

极小值点（1.0000，1.0000）
极小值y =1.9429e-16

求极值时也可以利用二阶导数：

编写目标函数fun4
function [f,df,d2f]=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)]; 
d2f=[-400*x(2)+1200*x(1)^2+2,-400*x(1) 
 -400*x(1),200]; 
编写主程序文件example22
options = optimset('GradObj','on','Hessian','on');
[x,y]=fminunc('fun4',rand(1,2),options)

• 例题3

% 求函数 f (x) = sin(x) + 3 取最小值时的 x 值
编写目标函数fun5
function f=fun5(x); 
f=sin(x)+3; 
编写主程序文件example3
x0=2; 
[x,y]=fminsearch(@fun5,x0) 

@fun5匿名函数详细介绍

2.2 求解方程的解和函数的零点

• 例题4求函数的零点

clc,clear
xishu=[1 -1 2 -3]; % 多项式是用向量定义的，系数从幂的高低排列
x0=roots(xishu) % 求得  -0.1378 ± 1.5273i和1.2757
 
%使用符号求解
syms x % 声明关键字
x0 = solve(x^3-x^2+2*x-3) % 求函数零点符号解
x0 = vpa(x0,5) % 化成小数格式
%  求得  -0.1378 ± 1.5273i和1.2757

%求数值解
y = @(x) x^3-x^2+2*x-3;
x = fsolve(y,rand) % 只能求解初始值附近的一个零点
% 求得   1.2757

• 例题5求方程的解

% 求下列方程组
x^2 + y - 6 = 0
y^2 + x - 6 = 0
% 符号求解方程解
syms s,y
[x,y] = solve(x^2 + y - 6,y^2 + x - 6) % 具体可以在Matlab中试一下输出值（一共四组解）

% 数值来求解方程的解
f = @(x)[x^2 + y - 6,y^2 + x - 6];
xy = fsolve(f,rand(2,1)) % 只能求解初始值附近的一个零点
% 求得一组数值解（2，2）

3.约束极值问题

3.1二次规划

若某非线性规划的目标函数为自变量x的二次函数，约束条件全是线性，就称为二次规划。
二次规划的数学模型：


其中H为实对称矩阵，f,b,beq,lb,ub为列向量；A,Aeq为相应维数的矩阵。
求解二次规划的命令为：

• 例题6
  

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)) 
% 求得
x1 = 1.9500
x2 = 1.0500
minf(x)=-11.0250

3.2罚函数法

罚函数法求解非线性规划问题的思想是，利用问题中的约束函数作出适当的罚函数，由此构造出带参数的增广目标函数，把问题转化为无约束非线性规划问题。

• 例题7
  

% 定义增广目标函数，编写函数test
function g=test(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);

% 利用Matlab的求矩阵的极大值和极小值函数编写test1
function g=test1(x); 
M=50000; 
f=x(1)^2+x(2)^2+8; 
g=f-M*sum(min([x';zeros(1,2)]))-M*min(x(1)^2-x(2),0)+M*abs(-x(1)-x(2)^2+2); 

% 也可以修改增广目标函数的定义，编写test2
function g=test2(x); 
M=50000; 
f=x(1)^2+x(2)^2+8; 
g=f-M*sum(min([x';zeros(1,2)]))-M*min(x(1)^2-x(2),0)+... 
 M*abs(-x(1)-x(2)^2+2); 

% 命令界面输入
[x,y]=fminunc('test几',rand(2,1)) 
% 由于限制只能求得一个局部最优解，使得每次运行结果都不一样

• Matlab中可以直接利用max、min和sum函数

3.3fminbnd函数

单变量非线性函数在区间上的极小值
min f(x),x属于[x1,x2]
Matlab的命令为：

[x,fval] = fminbnd(fun,x1,x2,options)
返回值为极小点x和函数的极小值，fun是用M文件定义的函数、匿名函数或Matlab中的单变量数学函数。

• 例题8

%求函数f(x) = (x-3)^2-1(x范围为[0,5])的最小值
编写函数example8
function f=example8(x)
f = (x-3)^2-1
% 命令界面输入
[x,y] = fminbnd('example8',0.5)
% 求得极小值点x=3,极小值为y=-1

3.4fseminf函数

Matlab的命令为：

[x,fval] = fseminf(fun,x0,ntheta,seminfcon,A,b,Aeq,beq,lb,ub)

ntheta是半无穷约束Ki(X,Wi)的个数；函数seminfcon用于定义非线性不等式c(x)、非线性等式约束ceq(x)和半无穷约束Ki(X,Wi)的函数，函数seminfcon有两个输入参数x和s。

3.5fminimax函数

• 例题9
  

编写函数example9
function f = example9(x)
f = [2*x(1)^2+x(2)^2-48*x(1)-40*x(2)+304
     -x(1)^2-3*x(2)^2
     x(1)+3*x(2)-18
     x(1)+x(2)-8];
% 调用函数fminimax
[x,y] = fminimax(@ example9,rand(2,1))

输出为

x =
    4.0000
    4.0000
y =
    0.0000
  -64.0000
   -2.0000
   -0.0000