matlab 向凸优化非线性约束函数传递参数 fmincon

您好，un为目标函数，它可用前面的方法定义；

x0为初始值；

A、b满足线性不等式约束 ，若没有不等式约束，则取A=[ ]，b=[ ]；

Aeq、beq满足等式约束 ，若没有，则取Aeq=[ ]，beq=[ ]；

lb、ub满足 ，若没有界，可设lb=[ ]，ub=[ ]；

nonlcon的作用是通过接受的向量x来计算非线性不等约束 和等式约束 分别在x处的估计C和Ceq，通过指定函数柄来使用，如：>>x = fmincon(@myfun,x0,A,b,Aeq,beq,lb,ub,@mycon)，先建立非线性约束函数，并保存为mycon.m：function [C,Ceq] = mycon(x)

C = …

% 计算x处的非线性不等约束 的函数值。

Ceq = …

% 计算x处的非线性等式约束 的函数值。

lambda是Lagrange乘子，它体现哪一个约束有效。

output输出优化猜陪首信息；

grad表示目标函数在x处的梯度；

hessian表示目标函数在x处的Hessiab值。

注意穗数：

1. fmincon 函数提供了大型优化算法和中型优化算法。默认时，若在 fun 函数中提供了梯度（options 参数的 GeadObj 设置为 'on'），并且只有上下界存在或只有等式约束，fmincon 函数将选择大型算法。 当既有等式约束又有梯度约束时，使用中型算法。

2. fmincon 函数的中型算法一般是使用序列二次规划。乱敏在每一步迭代中求解二次规划子问题，并用 BFGS 法更新 Lagrangian 乘子和 Hessian 矩阵。

3. fmincon 函数的大型算法采用了subspace trust region 优化算法。这种算法是把目标函数在点x的邻域泰勒展开（x可以认为是人为提供的初始猜测），这个展开的邻域就是所谓的trust region，泰勒展开进行到二阶项为止。

4. fmincon 函数可能会给出局部最优解，这与初始值的选取有关。

j是复数单位吗？可以试试http://dlib.net/，下面是个示例，我没有提供导数什么的，所以不稳定。它的用法和matlab挺像，你应该能搞定。

#include <dlib/optimization.h>

#include <iostream>

using namespace std

using namespace dlib

typedef matrix<double, 2, 1> column_vector

typedef complex<double> complex_t

int const M = 3

static matrix<double, 3, M> a

void init(){

    static bool tag = true

    if (tag){

        a = 0, 1, 2, 3, 4, 5, 6, 7, 8

        tag = false

    }

}

double foo(const column_vector& x){

    init()

    complex_t res

    for (int i 档谨= 0 i != a.nc() ++i){

        res += exp(complex_t(

            0,

            a(0, i) + a(1, i)*x(0) + a(2, i)*x(1)

            ))

    }

    return norm(res) / M

}

int main(){

    try{

        column_vector starting_point

        starting_point = 0, 0

        double max_value = find_max_box_constrained(bfgs_search_strategy(),

            objective_delta_stop_strategy(1e-9),

            foo, derivative(foo), starting_point, -50, 50)

        cout << "x0 = \n" << starting_point << "\nfoo(x0) = " << max_value <<氏肆 "\n"行核基

    }

    catch (std::exception const& e){

        cerr << e.what() << "\n"

    }

    return 0

}

运行结果：

用MATLAB优化工具箱解线性规划

命令：x=linprog（c，A，b）

命令：x=linprog（c，A，b，Aeq,beq）

注意：若没有不等式： 存在，则令A=[ ]，b=[ ]. 若没有等式约束, 则令Aeq=[ ], beq=[ ].

命令：[1] x=linprog（c，A，b，Aeq,beq, VLB，VUB）

[2] x=linprog（c，A，b，Aeq,beq, VLB，VUB, X0）

注意：[1] 若没有等式约束, 则令Aeq=[ ], beq=[ ]. [2]其中X0表示初始点

4、命令：[x,fval]=linprog(…)

返回最优解x及x处的目标函数值fval.

例1

解 编写M文件小xxgh1.m如下：

c=[-0.4 -0.28 -0.32 -0.72 -0.64 -0.6]

A=[0.01 0.01 0.01 0.03 0.03 0.030.02 0 0 0.05 0 00 0.02 0 0 0.05 00 0 0.03 0 0 0.08]

b=[850700100900]

Aeq=[]beq=[]

vlb=[000000]vub=[]

[x,fval]=linprog(c,A,b,Aeq,beq,vlb,vub)

例2 解: 编写M文件xxgh2.m如下：

c=[6 3 4]

A=[0 1 0]

b=[50]

Aeq=[1 1 1]

beq=[120]

vlb=[30,0,20]

vub=[]

[x,fval]=linprog(c,A,b,Aeq,beq,vlb,vub

例3 （任务分配问题）某车间有甲、乙两台机床，可用于加工三种工件。

假定这两台车床的可用台时数分别为800和900，三种工件的数量分别为400、

600和500，且已知用三种不同车床加工单位数量不同工件所需的台时数和加工

费用如下表。问怎样分配车床的加工任务，才能既满足加工工件的要求，又使

加工费用最低

解设在甲车床上加工工件1、2、3的数量分别为x1、x2、x3，在乙车床上

加工工件1、2、3的数量分别为x4、x5、x6。可建立以下线性规划模型：

编写M文件xxgh3.m如下:

f = [13 9 10 11 12 8]

A = [0.4 1.1 1 0 0 0

0 0 0 0.5 1.2 1.3]

b = [800900]

Aeq=[1 0 0 1 0 0

0 1 0 0 1 0

0 0 1 0 0 1]

beq=[400 600 500]

vlb = zeros(6,1)

vub=[]

[x,fval] = linprog(f,A,b,Aeq,beq,vlb,vub)

例4．某厂每日8小时的产量不低于1800件。为了进行质量控制，计划聘请两种不同水平的检验员。一级检验员的标准为册者：速度25件/小时，正确率98%，计时工资4元/小时；二级检验员的标准为：速度15小时/件，正确率95%，计时工资3元/小时。检验员每错检一次，工厂要损失2元。为使总检验费用最省，该工厂应聘一级、二级检验员各几名？

解设需要一级和二级检验员的人数分别为x1、x2人,

编写M文件xxgh4.m如下：

c = [4036]

A=[-5 -3]

b=[-45]

Aeq=[]

beq=[]

vlb = zeros(2,1)

vub=[915]

%调用linprog函数：

[x,fval] = linprog(c,A,b,Aeq,beq,vlb,vub)

结果为：

x =

9.0000

0.0000

fval =360

即只需聘用9个一级检验员。

4．控制参数options的设置

Options中常用的几个参数的名称、含义、取值如下:

(1) Display: 显示水平.取值为’off’时,不显示输出取值为’iter’时,显示每次迭代的信息取值为’final’时,显示最终结果.默认值为’final’.

(2) MaxFunEvals: 允粗塌许进行函数评价的最大次数,取值岩姿圆为正整数.

(3) MaxIter: 允许进行迭代的最大次数,取值为正整数

控制参数options可以通过函数optimset创建或修改。命令的格式如下：

(1) options=optimset(‘optimfun’)

创建一个含有所有参数名,并与优化函数optimfun相关的默认值的选项结构options.

（2）options=optimset(‘param1’,value1,’param2’,value2,...)

创建一个名称为options的优化选项参数,其中指定的参数具有指定值,所有未指定的参数取默认值.

(3)options=optimset(oldops,‘param1’,value1,’param2’,

value2,...)

创建名称为oldops的参数的拷贝,用指定的参数值修改oldops中相应的参数.

例：opts=optimset(‘Display’,’iter’,’TolFun’,1e-8)

该语句创建一个称为opts的优化选项结构,其中显示参数设为’iter’, TolFun参数设为1e-8.

用Matlab解无约束优化问题

一元函数无约束优化问题

常用格式如下：

（1）x= fminbnd (fun,x1,x2)

（2）x= fminbnd (fun,x1,x2 ，options)

（3）[x，fval]= fminbnd（...）

（4）[x，fval，exitflag]= fminbnd（...）

（5）[x，fval，exitflag，output]= fminbnd（...）

其中（3）、（4）、（5）的等式右边可选用（1）或（2）的等式右边。

函数fminbnd的算法基于黄金分割法和二次插值法，它要求目标函数必须是连续函数，并可能只给出局部最优解。

例1 求 在0<x<8中的最小值与最大值

主程序为wliti1.m:

f='2*exp(-x).*sin(x)'

fplot(f,[0,8])%作图语句

[xmin,ymin]=fminbnd (f, 0,8)

f1='-2*exp(-x).*sin(x)'

[xmax,ymax]=fminbnd (f1, 0,8)

运行结果：

xmin = 3.9270ymin = -0.0279

xmax = 0.7854 ymax = 0.6448

例2 对边长为3米的正方形铁板，在四个角剪去相等的正方形以制成方形无盖水槽，问如何剪法使水槽的容积最大？

先编写M文件fun0.m如下:

function f=fun0(x)

f=-(3-2*x).^2*x

主程序为wliti2.m:

[x,fval]=fminbnd('fun0',0,1.5)

xmax=x

fmax=-fval

运算结果为: xmax = 0.5000,fmax =2.0000.即剪掉的正方形的边长为0.5米时水槽的容积最大,最大容积为2立方米.

2、多元函数无约束优化问题

标准型为：min F(X)

命令格式为:

（1）x= fminunc（fun,X0 ）；或x=fminsearch（fun,X0 ）

（2）x= fminunc（fun,X0 ，options）；

或x=fminsearch（fun,X0 ，options）

（3）[x，fval]= fminunc（...）；

或[x，fval]= fminsearch（...）

（4）[x，fval，exitflag]= fminunc（...）；

或[x，fval，exitflag]= fminsearch

（5）[x，fval，exitflag，output]= fminunc（...）；

或[x，fval，exitflag，output]= fminsearch（...）

说明:

• fminsearch是用单纯形法寻优. fminunc的算法见以下几点说明：

[1] fminunc为无约束优化提供了大型优化和中型优化算法。由options中的参数LargeScale控制：

LargeScale=’on’(默认值),使用大型算法

LargeScale=’off’(默认值),使用中型算法

[2] fminunc为中型优化算法的搜索方向提供了4种算法，由

options中的参数HessUpdate控制：

HessUpdate=’bfgs’（默认值），拟牛顿法的BFGS公式；

HessUpdate=’dfp’，拟牛顿法的DFP公式；

HessUpdate=’steepdesc’，最速下降法

[3] fminunc为中型优化算法的步长一维搜索提供了两种算法，

由options中参数LineSearchType控制：

LineSearchType=’quadcubic’(缺省值)，混合的二次和三

次多项式插值；

LineSearchType=’cubicpoly’，三次多项式插

• 使用fminunc和 fminsearch可能会得到局部最优解.

例3 min f(x)=(4x12+2x22+4x1x2+2x2+1)*exp(x1)

1、编写M-文件 fun1.m:

function f = fun1 (x)

f = exp(x(1))*(4*x(1)^2+2*x(2)^2+4*x(1)*x(2)+2*x(2)+1)

2、输入M文件wliti3.m如下:

x0 = [-1, 1]

x=fminunc(‘fun1’,x0)

y=fun1(x)

3、运行结果:

x= 0.5000 -1.0000

y = 1.3029e-10

例4 Rosenbrock 函数 f（x1，x2）=100（x2-x12）2+(1-x1)2

的最优解（极小）为x*=（1，1），极小值为f*=0.试用

不同算法（搜索方向和步长搜索）求数值最优解.

初值选为x0=（-1.2 , 2）.

1. 为获得直观认识，先画出Rosenbrock 函数的三维图形,

输入以下命令：

[x,y]=meshgrid(-2:0.1:2,-1:0.1:3)

z=100*(y-x.^2).^2+(1-x).^2

mesh(x,y,z)

2. 画出Rosenbrock 函数的等高线图,输入命令：

contour(x,y,z,20)

hold on

plot(-1.2,2,' o ')

text(-1.2,2,'start point')

plot(1,1,'o')

text(1,1,'solution')

3.用fminsearch函数求解

输入命令:

f='100*(x(2)-x(1)^2)^2+(1-x(1))^2'

[x,fval,exitflag,output]=fminsearch(f, [-1.2 2])

运行结果:

x =1.00001.0000

fval =1.9151e-010

exitflag = 1

output =

iterations: 108

funcCount: 202

algorithm: 'Nelder-Mead simplex direct search'

4. 用fminunc 函数

(1)建立M-文件fun2.m

function f=fun2(x)

f=100*(x(2)-x(1)^2)^2+(1-x(1))^2

(2)主程序wliti44.m

Rosenbrock函数不同算法的计算结果

可以看出，最速下降法的结果最差.因为最速下降法特别不适合于从一狭长通道到达最优解的情况.

例5 产销量的最佳安排

某厂生产一种产品有甲、乙两个牌号，讨论在产销平衡的情况下如何确定各自的产量，使总利润最大. 所谓产销平衡指工厂的产量等于市场上的销量.

符号说明

z(x1,x2)表示总利润；

p1，q1，x1分别表示甲的价格、成本、销量；

p2，q2，x2分别表示乙的价格、成本、销量；

aij，bi，λi,ci（i，j =1，2）是待定系数.

基本假设

1．价格与销量成线性关系

利润既取决于销量和价格，也依赖于产量和成本。按照市场规律，

甲的价格p1会随其销量x1的增长而降低，同时乙的销量x2的增长也

会使甲的价格有稍微的下降，可以简单地假设价格与销量成线性关系，

即： p1 = b1 - a11 x1 - a12 x2 ，b1，a11，a12 >0，且a11 >a12；

同理， p2 = b2 - a21 x1- a22 x2 ，b2，a21，a22 >0

2．成本与产量成负指数关系

甲的成本随其产量的增长而降低,且有一个渐进值,可以假设为

负指数关系,

总利润为： z(x1,x2)=(p1-q1)x1+(p2-q2)x2

若根据大量的统计数据,求出系数b1=100,a11=1,a12=0.1,b2=280,

a21=0.2,a22=2,r1=30,λ1=0.015,c1=20, r2=100,λ2=0.02,c2=30,则

问题转化为无约束优化问题：求甲,乙两个牌号的产量x1，x2，使总利润z最大.

为简化模型,先忽略成本,并令a12=0,a21=0,问题转化为求:

z1 = ( b1 - a11x1 ) x1 + ( b2 - a22x2 ) x2

的极值. 显然其解为x1 = b1/2a11 = 50, x2 = b2/2a22 = 70,

我们把它作为原问题的初始值.

模型求解

1.建立M-文件fun.m:

function f = fun(x)

y1=((100-x(1)- 0.1*x(2))-(30*exp(-0.015*x(1))+20))*x(1)

y2=((280-0.2*x(1)- 2*x(2))-(100*exp(-0.02*x(2))+30))*x(2)

f=-y1-y2

2.输入命令:

x0=[50,70]

x=fminunc(‘fun’,x0),

z=fun(x)

3.计算结果:

x=23.9025, 62.4977, z=6.4135e+003

即甲的产量为23.9025,乙的产量为62.4977,最大利润为6413.5.

---