文章目录
 
1、LP问题何时局部最优解也是全局最优解
2、无约束非线性规划
2.1理论上的海塞矩阵、鞍点
2.2数值解法：最速下降法
    
*3、有约束非线性规划-KKT条件
*4、可行方向法
  
 

 
1、LP问题何时局部最优解也是全局最优解
 
对maximum问题来说：可行域是一个凸集，目标函数是凹规划
 对minimum问题来说：可行域是凸集，目标函数是凸规划
 
2、无约束非线性规划
 
2.1理论上的海塞矩阵、鞍点
 

 找局部最大值、局部最小值以及saddle point（鞍点），鞍点指的是不是极值点的驻点。一阶主子式应该是以-1的一次方为符号，二阶主子式应该是以-1的二次方为符号······这是concave凹函数，如果主子式的行列式的值全都是正的，那么就应该是convex凸函数。
 
 
2.2数值解法：最速下降法
 
 
*3、有约束非线性规划-KKT条件
 
使用kkt条件之前先判断限制条件的梯度之间是否是线性独立的。
 如果目标函数的梯度不可以使用条件约束的梯度表示，那么kkt条件不成立。不可以使用这种方法。
 
首先构建拉格朗日等式；
 然后分别对每个变量求解一阶导，令其等于0，求解出变量和拉朗日系数的值；
 判断目标函数和约束条件的凹凸性，然后得出局部最优解和全局最优解之间的关系；
 
*4、可行方向法

 
文章目录
 
数学建模笔记（三）非线性规划和多目标优化
非线性规划概念和理论
非线性规划问题的数学模型
有约束非线性规划的求解
凸规划
    
    
一个简单的非线性规划模型
二次规划模型
多目标规划问题
  
 

 
数学建模笔记（三）非线性规划和多目标优化
 
非线性规划概念和理论
 
非线性规划问题的数学模型
 
 
有约束非线性规划的求解
 
常见的处理思路是：可能的话将非线性问题转化为线性问题，将约束问题转化为无约束问题。
 1.对有等式约束的非线性规划问题使用Lagrange乘数法求解
 2.罚函数法：
 利用目标函数和约束函数构造带参数的增广目标函数，从而转换成一系列无约束非线性规划问题来进行求解。
 
凸规划
 
基本概念：
 1.凸集与凸函数的定义
 凸集：
 凸函数：
 例题1.
 求解：
 先对目标函数和约束函数中的非线性函数求二阶行列式，若均大于0，则是凸优化问题。
 求解代码：
 
import numpy as np
import cvxpy as cp
x=cp.Variable(2,pos=True)
obj=cp.Minimize(sum(x**2)-4*x[0]+4)
con=[-x[0]+x[1]-2<=0,
x[0]**2-x[1]+1<=0]
prob = cp.Problem(obj, con)
prob.solve(solver='CVXOPT')
print("最优值为:",round(prob.value,4))
print("最优解为：\n", np.round(x.value,4))
#-------------结果--------------
最优值为: 3.7989
最优解为：
 [0.5536 1.3064]
 
一个简单的非线性规划模型
 
数学建模五步骤：
 提出问题、选择建模方法、推导模型的数学表达式、求解模型、回答问题
 例题2.
 19英寸数量 $x_1$
 21英寸数量 $x_2$
 19英寸利润 $x_1\ast (339-195-0.01x_1-0.003x_2)$
 21英寸利润 $x_2\ast (399-225-0.01x_2-0.004x_1)$
 总利润 $x_1\ast (339-195-0.01x_1-0.003x_2)+x_2\ast (399-225-0.01x_2-0.004x_1)-400000$
 问题求解：
 

import sympy as sp
import pylab as plt
plt.rc('text', usetex=True) #使用 LaTeX字体
plt.rc('font',size=14)
sp.var('x1, x2') #定义符号变量
y = (339-0.01*x1-0.003*x2)*x1+(399-0.004*x1-0.01*x2)*x2-(400000+195*x1+225*x2)
y = sp.simplify(y) #化简
dy1 = y.diff(x1) #求关于 x1的偏导
dy2 = y.diff(x2) #求关于 x2的偏导
s = sp.solve([dy1, dy2], [x1, x2])
x10 = round(float(s[x1])) #取整
x20 = round(float(s[x2]))
y0 = y.subs({x1: x10, x2: x20}) #符号函数代入数值
f = sp.lambdify('x1, x2', y, 'numpy') #符号函数转换为匿名函数
x = plt.linspace(0, 10000, 100)
X, Y = plt.meshgrid(x, x) #转换为网格数据
Z = f(X, Y)
ax=plt.subplot(121, projection='3d') #第一个子窗口三维画图
ax.plot_surface(X, Y, Z,cmap='viridis')
ax.set_xlabel('$x_1$'); ax.set_ylabel('$x_2$')
plt.subplot(122) #激活第二个子窗口
contr=plt.contour(X,Y,Z,10) #10条等高线
plt.clabel(contr) #等高线标注
plt.ylabel('$x_2$',rotation=0)
plt.xlabel('$x_1$')
sp.var('a', pos=True) #定义灵敏度 分析 的符号参数
y = (339-a*x1-0.003*x2)*x1+(399-0.004*x1-0.01*x2)*x2-(400000+195*x1+225*x2)
y = sp.simplify(y) #化简
dy1 = y.diff(x1) #求关于 x1的偏导
dy2 = y.diff(x2) #求关于 x2的偏导
s = sp.solve([dy1, dy2], [x1, x2])
sx1 = s[x1]; sx2 = s[x2] #提取解分量
s1 = sp.lambdify('a', sx1, 'numpy') #符号函数转换为匿名函数
s2 = sp.lambdify('a', sx2, 'numpy')
a0 = plt.linspace(0.002, 0.02, 50)
plt.figure()
plt.subplots_adjust(wspace = 0.65)
plt.subplot(121); plt.plot(a0, s1(a0))
plt.xlabel('$a$'); plt.ylabel('$x_1$')
plt.subplot(122); plt.plot(a0, s2(a0))
plt.xlabel('$a$'); plt.ylabel('$x_2$')
#plt.subplot_tool() #调整位置的工具
dx1 = sx1.diff(a); dx10 = dx1.subs(a, 0.01)
sx1a = dx10 * 0.01 / 4735
dx2 = sx2.diff(a); dx20 = dx2.subs(a, 0.01)
sx2a = dx20 * 0.01 / 7043
Y = y.subs({x1: s[x1], x2: s[x2]}) #求关于 a的目标函数
Y = sp.factor(Y); Y = sp.simplify(Y)
Ya = sp.lambdify('a', Y, 'numpy') #转换为匿名函数
a0 = plt.linspace(0.002, 0.02, 1000)
plt.figure(); plt.plot(a0, Ya(a0))
plt.xlabel('$a$'); plt.ylabel('$y$', rotation=0)
Sya = - 4735 ** 2 * 0.01 / 553641.025
y2 = y.subs({x1: 4735, x2: 7043, a: 0.011}) #计算近似最优利润
y3 = Y.subs(a, 0.011) #计算最优利润
delta = (y3 - y2) / y2 #计算利润的相对误差
plt.show()
 
二次规划模型
 
定义：
 ### 非线性规划的求解和应用
 对于一般的非线性规划问题，由于不是凸优化，就不能使用cvxpy库求解。而需要使用scipy.optmize模块的minimize函数求解。
 例题5.
 由于不是凸函数，不能使用cvxpy库求解。
 求解程序：
 
import numpy as np
from scipy.optimize import minimize
c2 = np.array([[-1, -0.15],[-0.15, -2]])
c1 = np.array([98, 277])
a = np.array([[1, 1], [1, -2]])
b = np.array([100, 0])
obj = lambda x: x @ c2 @ x + c1 @ x
con ={'type': 'ineq', 'fun': lambda x: b-a@x}
bd = [(0, np.inf) for i in range(a.shape[1])]
res = minimize(obj, np.ones(2), constraints=con, bounds=bd)
print(res) #输出解的信息
#-------------
 fun: 0.0
     jac: array([ 97.99999999, 276.99999997])
 message: 'Optimization terminated successfully'
    nfev: 6
     nit: 2
    njev: 2
  status: 0
 success: True
       x: array([0., 0.])
 
多目标规划问题
 
定义：
 一般多目标规划问题没有最优解，只有满意解或有效解。
 求有效解的预处理：
 （1）无量纲化处理：每个目标函数的量纲通常是不一样的，在进行加权求解时由于量
 纲的不可公度性，需要先进行无量纲化处理。
 （2）数量级的归一化处理：当各个目标函数的数量级差异较大时，容易出现大数吃小
 数现象，即数量级较大的目标在决策分析过程中容易占优，从而影响决策结果。
 求解方法：
 1.线性加权法
 基本思想为根据目标的重要性确定一个权重，以目标函数的加权平均值为评价函数，是其达到最优。
 2.$ϵ$约束化
 其中${ϵ}_i$是决策者对第i个目标容许接受阈值。
 3.理想点法
 基本思想：以每个单目标最优值为该目标的理想值，使每个目标函数值与理
 想值的差的加权平方和最小。
 4.优先级法
 基本思想：根据目标重要性分成不同优先级，先求优先级高的目标函数的最优
 值，在确保优先级高的目标获得不低于最优值的条件下，再求优先级低的目标函数。

 
 
    
 


 

  

   
    与线性规划不同的是，非线性规划要求目标函数或约束条件中含有非线性函数。相应的求解这类问题就要用到非线性规划的方法。约束条件或者目标函数的放宽使得规划模型更具普适性，但也增加了问题求解的难度。对于简单的非线性规划问题，R语言中stat包即可求解。在这里我们给大家介绍R语言中求解非线性规划更为专业的Rdonlp2包。
  
 
 

  

   
  
 
 


 
    Rdonlp2在求解非线性规划问题上功能十分强大，用户可自行安装并查阅帮助文档：
 
library(Rdonlp2)

 
    Rdonlp2包的核心函数为 donlp2，可以快速求解非线性规划的最值。其用法如下：

 

 
donlp2(par,fn,
       par.upper=rep(+Inf,length(par)),
       par.lower=rep(+Inf,length(par)),
       A=NULL,
       lin.upper=rep(+Inf,length(par)),
       lin.lower=rep(-Inf,length(par)),
       nlin=list(),
       nlin.upper=rep(+Inf,length(nlin)),
       nlin.lower=rep(-Inf,length(nlin)),
       control=donlp2.control(),
       control.fun=function(lst){return(TRUE)},
       env=.GlobalEnv,
       name="Rdonlp2"
       )

 
    对该函数众多参数进行解释：

 
    par : 初始迭代向量

 
    fn : 连续型函数

 
    par.upper、par.lower : 自变量的上下限

 
    A：线性约束矩阵

 
    lin.upper、lin.lower : 线性约束条件的上下界

 
    nlin : 非线性约束条件函数列表

 
    nlin.upper、nlin.lower : 非线性约束条件的上下界

 
    其余参数可忽略。

 


 
    先看一例donlp2函数求解非线性规划的例子：   
 
    用Rdonlp2包编写非线性规划计算命令：
 
library(Rdonlp2)
p=c(0,0)              #迭代初始值
par.l=c(-100,100);par.u=c(100,100)#自变量定义域约束
fn=function(x){
   exp(x[1])*(4*x[1]^2+2*x[2]^2+4*x[1]*x[2]+2*x[2]+1)
}                  #目标函数
A=matrix(c(1,1,3,-1),2,byrow=TRUE)
lin.l=c(2,1);lin.u=c(+Inf,3)   #线性约束
nlcon1=function(x){
  -x[1]*x[2]
}                  #非线性约束
donlp2(p,fn,par.u,par.l,A,lin.l=lin.l,lin.u=lin.u,
       nlin=list(nlcon1))
       

$message
[1] "KT-conditions satisfied, no further correction computed"
$par
[1]  33.66667 100.00000
$gradf
[1]  1.625065e+19 2.243635e+17

 
    相应运行结果包括中间过程、算法参数、运行时间等数据，这里就不全部列出。

 
    再看一例利用Rdonlp2包求解二元rastrigin函数最小值的问题。该函数因为有众多极值点而仅有一个最值点对各类优化算法有一定的欺骗性，但用其来测试各类优化算法效果显著。该二元rastrigin函数为：
 

 
    先作出该函数的三维图像： 
 
a<-5
x<-seq(-a,a,0.01)
y<-seq(-a,a,0.01)
f<-function(x,y){
   x^2-10*sin(2*pi*x)+y^2-10*cos(2*pi*y)+20
}
z<-outer(x,y,f)
image(x,y,z,col=heat.colors(24))
library(rgl)
zorder<-rank(z)
persp3d(x,y,z,col=rainbow(as.integer(max(zorder)))[zorder])

 

 
fn<-function(x){
    f<-sum(x^2*cos(2*pi*x)+10)
}
par<-c(-500,500)
ret<-donlp2(par,fn)

ret$f
[1] 20
ret$par
[1] -2.654965e-06  2.654965e-06

 
    由输出结果可知该函数最小值为20 。所以Rdonlp2求解非线性规划也是非常方便的。当然了，还有很多更为一般的非线性规划问题是Rdonlp2包也不能解决的，这时候可能会向遗传算法、模拟退火这些启发式算法求助啦。这里且不作讨论。

 


 
    在许多实际问题中，衡量一个方案的好坏标准往往不止一个，例如制造一枚导弹，既要射程远又要最省燃料，还得精度高。这一类问题就不是单一的目标规划问题了，我们称之为多目标规划问题。

 
    求解多目标规划问题通常有主要目标法、分层序列法和线性加权求和法等方法。其中主要目标法通过确定一个主要目标，将多目标优化问题转化为线性或非线性规划问题；分层序列法是将目标中多个目标按照其重要程度排一个次序而逐步求解的过程；线性加权法则是对各目标赋予一个权数，加权求和得到一个新的目标函数而进行求解的过程。

 
    R语言中多目标规划问题的求解通常用goalprog包，其核心函数为llgp:
 
llgp(coefficient,targets,achievements,maxiter=1000,verbose=FALSE)

 
   其中coefficient为约束系数矩阵，targets为系数矩阵约束向量，achievements是包含objective、priority、p、n的目标函数，其中objective表示第几对偏差变量，priority表示该偏差变量的优先级，p和n分别为正负偏差变量的权系数。
 
    看一道多目标规划问题计算实例：（该题来自钱颂迪《运筹学》教材4.4节例5）

 

 
    llgp函数计算过程如下：
 
library(goalprog)
coefficients<-matrix(c(1,1,5,1,1,0,3,1),4)
targets<-c(10,4,56,12)
achievements<-data.frame(objective=1:4,priority=c(1,1,3,4),
                         p=c(2,3,0,1),n=c(0,0,1,0))
soln<-llgp(coefficients,targets,achievements)

soln$converged
[1] TRUE
soln$out
Decision variables
                X
X1   4.000000e+00
X2   6.000000e+00

Summary of objectives
        Objective           Over          Under         Target
G1   1.000000e+01   0.000000e+00   0.000000e+00   1.000000e+01
G2   4.000000e+00   0.000000e+00   0.000000e+00   4.000000e+00
G3   3.800000e+01   0.000000e+00   1.800000e+01   5.600000e+01
G4   1.000000e+01   0.000000e+00   2.000000e+00   1.200000e+01

Achievement function
                A
P1   0.000000e+00
P2   0.000000e+00
P3   1.800000e+01
P4   0.000000e+00

 
    由输出的结果我们可以看到x1、x2最优值分别为4和6。

 
    由该例我们可以看出，R语言求解多目标规划问题只需按照函数格式将规划的目标函数、约束条件套入格式即可，过程非常简单。

 


 
参考文献：
 
魏太云.R软件与最优化[C]//中国r语言会议.2008.
 


 


 


 


 


 

  

   

    

     
    
    

     
    
    

     扫码关注数据科学家养成记
    
    

   
  
  


 
 


 


 


 

 

 
下面是三个非线性规划领域的算法。课堂上给予了详细的讲解，在实践环节让学生编程实现，从而可以实验复杂一些的例子，加深对算法的理解。下面共有四个程序grad，simplelinesearch，bfgs和phr，全部使用MATLAB语言编写。这些代码远未完善，可修改余地很大，仅供教学之用。
 
function gradf=grad(hfun,x)
 
%GRAD  数值法求函数在给定点处的导数值(一元函数)或梯度(多元函数)
 
%   gradf = grad(hfun,x0)  hfun是函数句柄或内联函数，x0是一定点或一批点(按列)；返回
 
%   值gradf是函数在该点处的导数或梯度。
 
%   要求函数能对成批的点求函数值。比如：feval(hfun,X)返回一个与X同列的行向量，对应于以
 
%   X每一列作为函数自变量而求得的函数值。
 
%
 
%   Reference: 《最优化计算原理与算法程序设计》， 粟塔山等编著， 国防科技大学出版社
 
%  (湖南)， 2001.
 
%
 
%   $Author: WBC $    $Date: 2003/10/25 $
 
n = length(x);
 
h=1e-3;     % 数值法求梯度的步长
 
w1=zeros(n,1);%h/2
 
w2=w1; %-h/2
 
w3=w1; %h
 
w4=w1; %-h
 
for i=1:n
 
x(i)=x(i)+h/2;
 
w1(i)=hfun(x);
 
x(i)=x(i)-h;
 
w2(i)=hfun(x);
 
x(i)=x(i)-h/2;
 
w4(i)=hfun(x);
 
x(i)=x(i)+2*h;
 
w3(i)=hfun(x);
 
x(i)=x(i)-h;
 
end
 
gradf=(8*(w1-w2)-(w3-w4))/(6*h);
 
function [fv,x,lambda,exitflag]=simplelinesearch(hf,rho,l,u,lambda0,fv0,x0,g,d)
 
%LINESEARCH 简单线搜索
 
%输入参数：
 
%        hf -- 函数句柄，目标函数
 
%       rho -- 实标量，简单线搜索的参数，介于0到0.5之间的数
 
%         l -- 实标量，简单线搜索的参数，介于0到1之间的数
 
%         u -- 实标量，简单线搜索的参数，介于0到1之间的数，满足u>l
 
%    alpha0 -- 实标量，步长的初始值
 
%       fv0 -- 实标量，一维搜索的初始目标函数值，即 hf(x0+alpha_0*d)
 
%        x0 -- 实列向量，当前点
 
%         g -- 实列向量，函数hf在当前点处的梯度
 
%         d -- 实列向量，函数hf在当前点处的搜索方向
 
%输出参数：
 
%         fv --  实标量，一维搜索完成后的目标函数值，即 hf(x0+alpha*d)
 
%          x -- 实列向量，下一个点
 
%     lambda -- 实标量，可接受步长
 
%   exitflag -- 整型标量，等于0表示线搜索成功，等于-1表示线搜索失败(内部迭代次数大于iterMax)
 
%参考文献：倪勤，最优化方法与程序设计，科学出版社
 
% Date: 2009/12/20
 
lambda = lambda0;
 
x = x0;
 
i = 0;
 
imax = 30; % 最大迭代次数，用户可以修改
 
% 主循环
 
gd = dot(g, d);
 
while i <= imax
 
fv = hf(x + lambda*d);
 
i = i + 1;
 
if fv < fv0 + lambda * rho * gd;
 
x = x + lambda*d;
 
exitflag = 0;
 
return;
 
end
 
lambda_bar = -gd*lambda^2*0.5/(fv-fv0-lambda*gd);
 
lambda_bar = max(lambda_bar, l*lambda);
 
lambda = min(lambda_bar, u*lambda);
 
end
 
exitflag=0;
 
if i >= imax && fv >= fv0
 
fv = fv0;
 
x = x0;
 
lambda = 0;
 
exitflag = -1;
 
end
 
function [fv, x, exitflag] = bfgs(hf, x0, epsi)
 
%BFGS 无约束问题的BFGS算法
 
%输入参数：
 
%        hf -- 函数句柄，目标函数
 
%        x0 -- 实列向量，初始点
 
%      epsi -- 实标量，终止误差
 
%输出参数：
 
%         fv --  实标量，一维搜索完成后的目标函数值，即 hf(x0+alpha*d)
 
%          x -- 实列向量，下一个点
 
%   exitflag -- 整型标量，等于0表示成功，等于-1表示失败(迭代次数大于iter_max)
 
%参考文献：倪勤，最优化方法与程序设计，科学出版社
 
% Date: 2009/12/20
 
%
 
%初始化
 
%
 
k = 0;
 
rho = 0.01;
 
l = 0.15;
 
u = 0.85;
 
x =x0;
 
fv = hf(x);
 
n = length(x);
 
H = eye(n);
 
iter_max = 100;
 
%
 
%检查终止条件
 
%
 
g = grad(hf, x);
 
while norm(g) > epsi && k<= iter_max
 
d = -H*g;
 
lambda0 = 1.0;
 
%
 
%做线搜索，如果成功，则返回更新当前点
 
%
 
[fv, x, lambda, exitflag] = simplelinesearch(hf, rho, l, u, lambda0, fv, x, g, d);
 
if exitflag == -1 %重开始
 
H = eye(n);
 
else
 
%
 
%更新H
 
%
 
g_old = g;
 
g = grad(hf, x);
 
p = lambda*d;
 
q = g- g_old;
 
Hq = H*q;
 
pq = p'*q;
 
qHq = q'*Hq;
 
v = sqrt(qHq) * (p/pq-Hq/qHq);
 
H = H + p*p'/pq - Hq*Hq'/qHq + v*v';
 
end
 
k = k+1;
 
end
 
exitflag = 0;
 
if k > iter_max
 
exitflag = -1;
 
end
 
function [fv, x, exitflag] = phr(hf, cf, x0)
 
%PHR  一般约束优化问题的PHR算法
 
%输入参数：
 
%         hf -- 函数句柄，目标函数
 
%         cf -- 函数句柄，约束条件，包含等式约束和不等式约束
 
%         x0 -- 实列向量，初始点
 
%输出参数：
 
%         fv --  实标量，一维搜索完成后的目标函数值，即 hf(x0+alpha*d)
 
%          x -- 实列向量，下一个点
 
%   exitflag -- 整型标量，等于0表示成功，等于-1表示失败(迭代次数大于iter_max)
 
%参考文献：倪勤，最优化方法与程序设计，科学出版社
 
% Date: 2009/12/20
 
%
 
%初始化
 
%
 
epsi = 1.0e-4;
 
k = 0;
 
sigma = 0.8;
 
c = 1.5;
 
theta = 0.8;
 
x = x0;
 
[ce, ci] = cf(x);
 
l = length(ce);
 
li = length(ci);
 
lambda = ones(l+li, 1) * 0.1;
 
iter_max = 100;
 
phi = 0;
 
if l
 
phi = phi + ce'*ce;
 
end
 
if li
 
phi = phi + sum(min(ci,lambda(l+1:end)/sigma).^2);
 
end
 
while phi > epsi && k <= iter_max
 
%
 
%求解子问题
 
%
 
hmf = @(x) mfun(x, hf, cf,lambda, sigma);
 
[fv, x] = bfgs(hmf, x, epsi);
 
[ce, ci] = cf(x);
 
phi_old = phi;
 
phi = 0;
 
if l
 
phi = phi + ce'*ce;
 
end
 
if li
 
phi = phi + sum(min(ci,lambda(l+1:end)/sigma).^2);
 
end
 
if phi > epsi
 
%
 
%更新罚因子
 
%
 
if k >= 2 && phi/phi_old > theta
 
sigma = c * sigma;
 
end
 
%
 
%更新乘子
 
%
 
if l
 
lambda(1:l) = lambda(1:l) - sigma*ce;
 
end
 
if li
 
lambda(l+1:end) = max(0, lambda(l+1:end) - sigma*ci);
 
end
 
end
 
k = k+1;
 
end
 
exitflag = 0;
 
fv = hf(x);
 
if k > iter_max
 
exitflag = -1;
 
end
 
%
 
%乘子罚函数
 
%
 
function fv = mfun(x, hf, cf, lambda, sigma)
 
[ce, ci] = cf(x);
 
l = length(ce);
 
li = length(ci);
 
fv = 0;
 
fv = fv + hf(x);
 
if l
 
fv = fv - lambda(1:l)'*ce + 0.5*sigma*ce'*ce;
 
end
 
if li
 
fv = fv + 0.5/sigma*sum(max(0,lambda(l+1:end) - sigma*ci).^2 - lambda(l+1:end).^2);
 
end
 
这里是演示代码：
 
%% bfgs演示
 
%教材P328.1-3
 
hf1 = @(x) 100 * (x(2) - x(1).^2).^2 + (1 - x(1)).^2; %banana函数
 
hf2 = @(x) (6 + x(1) + x(2)).^2 + (2 - 3*x(1) - 3*x(2) - x(1)*x(2)).^2;
 
hf3 = @(x) x(1).^2 - 2*x(1)*x(2) + 4*x(2).^2 + x(1) - 3*x(2);
 
[fv,x,exitflag]=bfgs(hf1,[0;0],0.001);
 
fv
 
x
 
exitflag
 
[fv,x,exitflag]=bfgs(hf2,[4;6],0.001);
 
fv
 
x
 
exitflag
 
[fv,x,exitflag]=bfgs(hf3,[1;1],0.001);
 
fv
 
x
 
exitflag
 
%% phr算法
 
%教材P414.5
 
hf1 = @(x) x(1).^2 + x(2).^2;
 
cf1 = @(x) deal([], x(1) - 1);
 
hf2 = @(x) x(1) + (x(2) + 1).^2/3;
 
cf2 = @(x) deal([],[x(1); x(2) - 1]);
 
%教材P392.2
 
hf3 = @(x) x(1).^2 + x(1).*x(2) + 2*x(2).^2 - 6*x(1) - 2*x(2) - 12*x(3);
 
cf3 = @(x) deal(x(1)+x(2)+x(3)-2,...
 
[x(1) - 2*x(2) + 3; x(1); x(2); x(3)]);
 
%其它例子
 
hf4 = @(x) 6*x(2)*x(5) + 7*x(1)*x(3) + 3*x(2)^2;
 
cf4 = @(x) deal([3*x(2)^2*x(5) + 3*x(1)^2*x(3) - 20.875;
 
x(1) - 0.3*x(2)],...
 
[-x(1) + 0.2*x(2)*x(5) + 71
 
-0.9*x(3) + x(4)^2 + 67
 
x(3)
 
x(5) - 1
 
-x(3) + 20
 
x(4) - 0.1*x(5)
 
-x(4) + 0.5*x(5)
 
x(3) - 0.9*x(5)
 
]);
 
hf5 = @(x) exp(x(1)) * (4*x(1)^2 + 2*x(2)^2 + 4*x(1)*x(2) + 2*x(2) + 1);
 
cf5 = @(x) deal([], [x(1) + x(2) - x(1)*x(2) - 1.5; x(1)*x(2) + 10]);
 
[fv, x, exitflag] = phr(hf1, cf1, [3;2]); fv x exitflag [fv, x, exitflag] = phr(hf2, cf2, [3;2]); fv x exitflag [fv, x, exitflag] = phr(hf3, cf3, [1;1;0]); fv x exitflag [fv, x, exitflag] = phr(hf4, cf4, [1; 4; 5; 2; 5]); fv x exitflag [fv, x, exitflag] = phr(hf5, cf5, [-1; 1]); fv x exitflag