机器学习深度学习高阶内容系列-Python实现凸优化求解器
 
CVXOPT是一个基于Python编程语言的凸优化的免费软件包。它可以与交互式Python解释器一起使用，通过执行Python脚本在命令行上使用，或者通过Python扩展模块集成到其他软件中。它的主要目的是通过构建Python广泛的标准库和Python作为高级编程语言的优势，使凸优化应用程序的软件开发变得简单明了。
 
CVXOPT用两个矩阵对象扩展了内置Python对象：一个用于密集矩阵的矩阵对象和一个用于稀疏矩阵的spmatrix对象。密集矩阵和稀疏矩阵一章介绍了这两种矩阵类型，以及为它们定义的算术运算和函数。 BLAS接口和稀疏线性方程 描述了几个用于密集和稀疏矩阵计算的库的接口。CVXOPT优化例程在Cone编程和建模章节中进行了描述。这些工具包括用Python编写的凸优化求解器、与其他一些优化库的接口，以及用于分段线性凸优化问题的建模工具。
 
 
 
 
目录 
 
 

  
CVXOPT 安装部署
  
创建矩阵
  
矩阵索引
  
Numpy和CVXOPT
  
求解线性规划
  
求解二次规划
 
 
 
 
CVXOPT 安装部署
 
官网链接：http://cvxopt.org/userguide/index.html
 

1 简介
 
函数优化是遗传算法的经典应用领域,也是遗传算法进行性能评价的常用算例.对于一些非线性,多模型,多目标的函数优化问题,用其它优化方法较难求解,而用遗传算法则可以方便地得到较好的结果.用Matlab实现了函数优化的遗传算法源程序.
 
遗传算法的流程如图 1 所示。其主要组成部分包括：　　
 
​（1）编码方式。遗传算法通常根据问题本身进行编码，并将问题的有效解决方案转化为遗传算法的搜索空间。工业中常用的编码方法包括实数编码，二进制编码，整数编码和数据结构编码。　
 
（2）适应度函数。适应度函数，也称为目标函数，是对整个个体与其适应度之间的对应关系的描述。具有高适应性的个体中包含的高质量基因具有较高的传递给后代的概率，而具有低适应性的个体的遗传概率较低。　　
 
（3）遗传操作。基本的遗传操作包括：选择、交叉、变异。　　
 
a）选择。选择操作基于个体适应度评估，选择群体中具有较高适应度的个体，并且消除具有较低适应度的个体。当然不同的选择操作也会带来不同的结果，有效的选择操作可以显著的提高搜索的效率和速度，减少无用的计算量。　　
 
常见的选择方法有：基于比例的适应度分配方法，期望值选择方法，基于排名的适应度分配方法，轮盘赌选择方法等。　　
 
b）交叉。在自然界生物进化过程中，两条染色体通过基因重组形成新的染色体，因此交叉操作是遗传算法的核心环节。交叉算子的设计需要根据具体的问题具体分析，编码操作和交叉操作互相辅助，交叉产生的新的个体必须满足染色体的编码规律。父代染色体的优良性状最大程度上的遗传给下一代染色体，在此期间也能能够产生一些较好的性状。　　
 
常见的交叉算子包括实质重组，中间重组，离散重组，线性重组，二进制交叉，单点交叉，均匀交叉，多点交叉和减少代理交叉。　　
 
c) 变异。通过随机选择的方法改变染色体上的遗传基因。变异本身可以被视为随机算法，严格来说，是用于生成新个体的辅助算法。　　
 
几个与浮点数编码和二进制编码个体匹配的交叉运算：单点交叉，均匀交叉，算术交叉，两点交叉和多点交叉。
 
 
2 部分代码
 
function varargout = yichuan(varargin)
% YICHUAN MATLAB code for yichuan.fig
%     YICHUAN, by itself, creates a new YICHUAN or raises the existing
%     singleton*.
%
%     H = YICHUAN returns the handle to a new YICHUAN or the handle to
%     the existing singleton*.
%
%     YICHUAN('CALLBACK',hObject,eventData,handles,...) calls the local
%     function named CALLBACK in YICHUAN.M with the given input arguments.
%
%     YICHUAN('Property','Value',...) creates a new YICHUAN or raises the
%     existing singleton*. Starting from the left, property value pairs are
%     applied to the GUI before yichuan_OpeningFcn gets called. An
%     unrecognized property name or invalid value makes property application
%     stop. All inputs are passed to yichuan_OpeningFcn via varargin.
%
%     *See GUI Options on GUIDE's Tools menu. Choose "GUI allows only one
%     instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help yichuan

% Last Modified by GUIDE v2.5 15-Jul-2016 21:52:11

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                  'gui_Singleton',  gui_Singleton, ...
                  'gui_OpeningFcn', @yichuan_OpeningFcn, ...
                  'gui_OutputFcn',  @yichuan_OutputFcn, ...
                  'gui_LayoutFcn', [] , ...
                  'gui_Callback',   []);
if nargin && ischar(varargin{1})
   gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
  [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
   gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT


% --- Executes just before yichuan is made visible.
function yichuan_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject   handle to figure
% eventdata reserved - to be defined in a future version of MATLAB
% handles   structure with handles and user data (see GUIDATA)
% varargin   command line arguments to yichuan (see VARARGIN)

% Choose default command line output for yichuan
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes yichuan wait for user response (see UIRESUME)
% uiwait(handles.figure1);


% --- Outputs from this function are returned to the command line.
function varargout = yichuan_OutputFcn(hObject, eventdata, handles) 
% varargout cell array for returning output args (see VARARGOUT);
% hObject   handle to figure
% eventdata reserved - to be defined in a future version of MATLAB
% handles   structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;


% --- Executes on button press in pushbutton1.
function pushbutton1_Callback(hObject, eventdata, handles)
% hObject   handle to pushbutton1 (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles   structure with handles and user data (see GUIDATA)
global x
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%绘制题中函数图形%%%%%%%%%%%%%%%%%%%%%%%%%%%
x=linspace(-1,2);
y=10+x.*cos(5*pi*x);

axes(handles.axes1)
plot(x,y,'r')
title('函数曲线图')
xlabel('x')
ylabel('y')


% --- Executes on button press in pushbutton2.
function pushbutton2_Callback(hObject, eventdata, handles)
% hObject   handle to pushbutton2 (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles   structure with handles and user data (see GUIDATA)
global BitLength         %全局变量，计算如果满足求解精度至少需要编码的长度 
global boundsbegin       %全局变量，自变量的起始点 
global boundsend         %全局变量，自变量的终止点 
=
end
Generation=Generation-1;%Generation加1、减1的操作是为了能记录各代中的最佳函数值xmax(Generation)
targetfunvalue=targetfun(xmax)  
[Besttargetfunvalue,nmax]=max(targetfunvalue) 
Bestpopulation=xmax(nmax)   
%绘制经过遗传运算后的适应度曲线
axes(handles.axes2)
hand1=plot(1:Generation,ymax);
set(hand1,'linestyle','-','linewidth',1,'marker','*','markersize',8) 
hold on;  

hand2=plot(1:Generation,ymean); 
set(hand2,'color','k','linestyle','-','linewidth',1, 'marker','h','markersize',8) 
xlabel('进化代数');
ylabel('最大和平均适应度');
xlim([1 Generationnmax]);
legend('最大适应度','平均适应度');
box off;
hold off;    


% --- Executes on button press in pushbutton3.
function pushbutton3_Callback(hObject, eventdata, handles)
% hObject   handle to pushbutton3 (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles   structure with handles and user data (see GUIDATA)
clearall;  
close all;
 
3 仿真结果
 
 
4 参考文献
 
[1]宋建萍. "函数优化的遗传算法代码实现." 软件导刊 12.002(2013):40-42.
 
部分理论引用网络文献，若有侵权联系博主删除。
 
5 MATLAB代码与数据下载地址
 
见博客主页

【优化求解】基于粒子群算法求解多目标优化问题matlab源码
 
 
 
1 算法介绍
 
1.1 关于速度和位置
 
粒子群算法通过设计一种无质量的粒子来模拟鸟群中的鸟，粒子仅具有两个属性：速度和位置，速度代表移动的快慢，位置代表移动的方向。
 
鸟被抽象为没有质量和体积的微粒(点)，并延伸到N维空间，粒子i在N维空间的位置表示为矢量Xi＝(x1，x2，…，xN)，飞行速度表示为矢量Vi＝(v1，v2，…，vN)。每个粒子都有一个由目标函数决定的适应值(fitness value)，并且知道自己到目前为止发现的最好位置(pbest)和现在的位置Xi。这个可以看作是粒子自己的飞行经验。除此之外，每个粒子还知道到目前为止整个群体中所有粒子发现的最好位置(gbest)(gbest是pbest中的最好值)，这个可以看作是粒子同伴的经验。粒子就是通过自己的经验和同伴中最好的经验来决定下一步的运动。
 
2.2 速度和位置的更新
 
PSO初始化为一群随机粒子(随机解)。然后通过迭代找到最优解。在每一次的迭代中，粒子通过跟踪两个“极值”(pbest，gbest)来更新自己。在找到这两个最优值后，粒子通过下面的公式来更新自己的速度和位置。
 
 
对于公式(1)：
 
公式(1)的第①部分称为【记忆项】，表示上次速度大小和方向的影响；
 
公式(1)的第②部分称为【自身认知项】，是从当前点指向粒子自身最好点的一个矢量，表示粒子的动作来源于自己经验的部分；
 
公式(1)的第③部分称为【群体认知项】，是一个从当前点指向种群最好点的矢量，反映了粒子间的协同合作和知识共享。粒子就是通过自己的经验和同伴中最好的经验来决定下一步的运动。
 
以上面两个公式为基础，再来看一个公式：
 
 
公式(2)和 公式(3)被视为标准PSO算法。
 
1.3 标准PSO算法的流程
 
1）初始化一群微粒(群体规模为N)，包括随机位置和速度；
 
2）评价每个微粒的适应度；
 
3）对每个微粒，将其适应值与其经过的最好位置pbest作比较，如果较好，则将其作为当前的最好位置pbest；
 
4）对每个微粒，将其适应值与其经过的最好位置gbest作比较，如果较好，则将其作为当前的最好位置gbest；
 
5）根据公式(2)、(3)调整微粒速度和位置；
 
6）未达到结束条件则转第2）步。
 
迭代终止条件根据具体问题一般选为最大迭代次数Gk或(和)微粒群迄今为止搜索到的最优位置满足预定最小适应阈值。
 
## 4.2 PSO流程图解
 
 
2 部分代码
 
 
clc;
clear;
close all;
​
%% 问题定义，这里可以删除
​
TestProblem=1;   % Set to 1, 2, or 3
​
switch TestProblem
    case 1
        CostFunction=@(x) MyCost1(x);
        nVar=50;
        VarMin=0;
        VarMax=1;
        
    case 2
        CostFunction=@(x) MyCost2(x);
        nVar=3;
        VarMin=-5;
        VarMax=5;
        
    case 3
        CostFunction=@(x) MyCost3(x);
        nVar=2;
        VarMin=0;
        VarMax=1;
end
​
VarSize=[1 nVar];
​
VelMax=(VarMax-VarMin)/10;
​
%% MOPSO 设置
nPop=500;   % Population Size
​
nRep=100;   % Repository Size
​
MaxIt=500;  % Maximum Number of Iterations
​
phi1=2.05;%学习因子
phi2=2.05;
phi=phi1+phi2;
chi=2/(phi-2+sqrt(phi^2-4*phi));
​
w=chi;              % Inertia Weight
wdamp=1;            % Inertia Weight Damping Ratio
c1=chi*phi1;        % Personal Learning Coefficient
c2=chi*phi2;        % Global Learning Coefficient
​
alpha=0.1;  % Grid Inflation Parameter
​
nGrid=10;   % Number of Grids per each Dimension
​
beta=4;     % Leader Selection Pressure Parameter
​
gamma=2;    % Extra (to be deleted) Repository Member Selection Pressure
​
%% 初始化
​
particle=CreateEmptyParticle(nPop);
​
for i=1:nPop
    particle(i).Velocity=0;
    particle(i).Position=unifrnd(VarMin,VarMax,VarSize); %unifrnd在[VarMin,VarMax]区间模拟VarSize的数组
    % 初始化了自变量X的位置
    particle(i).Cost=CostFunction(particle(i).Position); % 求解目标函数
    particle(i).Best.Position=particle(i).Position;
    particle(i).Best.Cost=particle(i).Cost;
end
​
particle=DetermineDomination(particle);
​
rep=GetNonDominatedParticles(particle);
​
rep_costs=GetCosts(rep);
G=CreateHypercubes(rep_costs,nGrid,alpha);
​
for i=1:numel(rep)
    [rep(i).GridIndex rep(i).GridSubIndex]=GetGridIndex(rep(i),G);
end
    
%% MOPSO 主循环
​
for it=1:MaxIt
    for i=1:nPop
        rep_h=SelectLeader(rep,beta);
​
        particle(i).Velocity=w*particle(i).Velocity ...
                             +c1*rand*(particle(i).Best.Position - particle(i).Position) ...
                             +c2*rand*(rep_h.Position -  particle(i).Position);
​
        particle(i).Velocity=min(max(particle(i).Velocity,-VelMax),+VelMax);
​
        particle(i).Position=particle(i).Position + particle(i).Velocity;
    
        flag=(particle(i).Position<VarMin | particle(i).Position>VarMax);    %非劣检查:flag
        particle(i).Velocity(flag)=-particle(i).Velocity(flag);
        
        particle(i).Position=min(max(particle(i).Position,VarMin),VarMax);
​
        particle(i).Cost=CostFunction(particle(i).Position);
​
        if Dominates(particle(i),particle(i).Best)
            particle(i).Best.Position=particle(i).Position;
            particle(i).Best.Cost=particle(i).Cost;
            
        elseif ~Dominates(particle(i).Best,particle(i))
            if rand<0.5
                particle(i).Best.Position=particle(i).Position; 
                particle(i).Best.Cost=particle(i).Cost;
            end
        end
​
    end
    
    particle=DetermineDomination(particle);
    nd_particle=GetNonDominatedParticles(particle);
    
    rep=[rep
         nd_particle];
    
    rep=DetermineDomination(rep);
    rep=GetNonDominatedParticles(rep);
    
    for i=1:numel(rep)
        [rep(i).GridIndex rep(i).GridSubIndex]=GetGridIndex(rep(i),G);
    end
    
    if numel(rep)>nRep
        EXTRA=numel(rep)-nRep;
        rep=DeleteFromRep(rep,EXTRA,gamma);
        
        rep_costs=GetCosts(rep);
        G=CreateHypercubes(rep_costs,nGrid,alpha);
        
    end
   
%     disp(['Iteration ' num2str(it) ': Number of Repository Particles = ' num2str(numel(rep))]);
    
    w=w*wdamp;
end
​
%% 结果
​
costs=GetCosts(particle);
rep_costs=GetCosts(rep);
​
figure;
​
plot(costs(1,:),costs(2,:),'b.');
hold on;
plot(rep_costs(1,:),rep_costs(2,:),'rx');
legend('Main Population','Repository');
​
 
 
3 仿真结果
 
 
 
4 参考文献
 
[1]王越, 吕光宏. 改进的粒子群求解多目标优化算法[J]. 计算机技术与发展, 2014, 000(002):42-45.
 
5 代码下载
 
 
 

优化工具包—无约束非线性优化求解器（fminunc）
 
原创不易，路过的各位大佬请点个赞
 


 
目录
 
优化工具包—无约束非线性优化求解器（fminunc）
一、fminunc总体介绍
二、fminunc求解器的具体用法
三、举例：最小化多项式
四、举例—获取最佳目标函数值
五、检查求解过程（options 设置）
六、信赖域法实例（fminunc求解器trust-region设置）
6.1 信赖域实例1
6.2 信赖域实例2：使用问题结构体
  
 


 
fminunc函数：求无约束多变量函数的最小值
 
 
一、fminunc总体介绍
 
语法：
 
x = fminunc(fun,x0)
x = fminunc(fun,x0,options)
x = fminunc(problem)
[x,fval] = fminunc(___)
[x,fval,exitflag,output] = fminunc(___)
[x,fval,exitflag,output,grad,hessian] = fminunc(___)
 
求解如下目标函数的最小值
 $$\underset{x}{\min }f(x)$$
 其中$f(x)$目标损失函数，即优化函数。$x$为变量，fminunc求解器面向多变量函数，因此$x$可以是向量和矩阵。与fminbnd求解器不同，它的变量必须是单变量。此外，fminunc 适用于无约束非线性问题。如果需要优化的问题有约束，通常使用 fmincon。
 
 
 
x = fminunc(fun,x0) 在点 x0 处开始并尝试求 fun 中描述的函数的局部最小值 x。点 x0 可以是标量、向量或矩阵。
 
 
二、fminunc求解器的具体用法
 
下面具体介绍fminunc求解器的用法
 
 
 
1）x = fminunc(fun,x0,options) 使用 options 中指定的优化选项最小化 fun。使用 optimoptions 可设置这些选项（后面介绍options）;
 
 
 
 
2）x = fminunc(problem) 求 problem 的最小值，它是 problem 中所述的一个结构体，描述一个优化问题；
 
 
 
 
3）[x,fval] = fminunc(___) 对上述任何语法，返回目标函数 fun 在解 x 处的值；
 
 
 
 
4）[x,fval,exitflag,output] = fminunc(___) 还返回描述 fminunc 的退出条件的值 exitflag，以及提供优化过程信息的结构体 output。
 
 
 
 
 
5）[x,fval,exitflag,output,grad,hessian] = fminunc(___) 还返回：
 ~grad - fun 在解 x 处的梯度。
 ~hessian - fun 在解 x 处的 Hessian 矩阵。请参阅 fminunc Hessian 矩阵。
 
 
 
三、举例：最小化多项式
 
目标函数：
 $$f(x)=3x_1^2+2x_2{x}_2+x_2^2-4x_1+5x_2$$
 其中$x=[x_1,x_2]$为一个二维向量。
 代码：
 给定初始值x0=[1,1];
 
fun = @(x)3*x(1)^2 + 2*x(1)*x(2) + x(2)^2 - 4*x(1) + 5*x(2);
x0 = [1,1];
[x,fval] = fminunc(fun,x0)
 
求解结果：
 
 目标函数的三维网格图：
 
fun = @(x) 3*x(1)^2 + 2*x(1)*x(2) + x(2)^2 - 4*x(1) + 5*x(2);
x0 = [1,1];
[x,fval] = fminunc(fun,x0)
xmin=x(1);ymin=x(2);zmin=fval;

syms x y
z=3*x^2 + 2*x*y + y^2 - 4*x + 5*y;
ezsurf(z,[-20,20],30);hold on;
plot3(xmin,ymin,zmin,'rp','markersize',30,'MarkerFaceColor','r')   %标记一个黑色的圆点
text(xmin,ymin,zmin,[' xmax=',num2str(xmin),char(10),'  y=',num2str(ymin),char(10),'  z=',num2str(zmin)]);   %标出坐标
 
 
当然也可以输出x处的梯度信息:
 
fun = @(x) 3*x(1)^2 + 2*x(1)*x(2) + x(2)^2 - 4*x(1) + 5*x(2);
x0 = [1,1];

[x,fval,exitflag,output,grad,hessian] = fminunc(fun,x0)
 
 
四、举例—获取最佳目标函数值
 
求非线性函数最小值的位置和该最小值处的函数值。目标函数是
 
 
fun = @(x)x(1)*exp(-(x(1)^2 + x(2)^2)) + (x(1)^2 + x(2)^2)/20;
x0 = [1,2];
[x,fval,exitflag,output,grad,hessian] = fminunc(fun,x0)
 
结果如下
 
 
 
x = -0.6691 0.0000
 fval = -0.4052
 exitflag =
 1 output = 包含以下字段的 struct:
 iterations: 9
 funcCount: 42
 stepsize: 2.9345e-04
 lssteplength: 1
 firstorderopt: 8.0839e-07
 algorithm: ‘quasi-newton’
 message: ‘Local minimum found.↵↵Optimization completed because the size of the gradient is less than↵the default value of the
 optimality tolerance.↵↵Stopping criteria details:↵↵Optimization
 completed: The first-order optimality measure, 6.891345e-07, is less
 ↵than options.OptimalityTolerance = 1.000000e-06.↵↵Optimization Metric
 Options↵relative norm(gradient) = 6.89e-07
 OptimalityTolerance = 1e-06 (default)’
 grad =
 1.0e-06 *
 0.8084
 0.6817
 hessian =
 1.8998 -0.0000 -0.0000 0.9552
 三维网格图如下：
 
 
close all;
fun = @(x)x(1)*exp(-(x(1)^2 + x(2)^2)) + (x(1)^2 + x(2)^2)/20;
x0 = [1,2];
[x,fval,exitflag,output,grad,hessian] = fminunc(fun,x0)

figure
xmin=x(1);ymin=x(2);zmin=fval;
syms x y
z=x*exp(-(x^2 + y^2)) + (x^2 + y^2)/20;
ezsurf(z,[-20,20],30);hold on;
plot3(xmin,ymin,zmin,'rp','markersize',30,'MarkerFaceColor','r')   %标记一个黑色的圆点
text(xmin,ymin,zmin,[' xmax=',num2str(xmin),char(10),'  y=',num2str(ymin),char(10),'  z=',num2str(zmin)]);   %标出坐标
 

 画出等高线
 
close all
figure
X=-20:0.1:20;
[x,y]=meshgrid(X);
z=x.*exp(-(x.^2 + y.^2)) + (x.^2 + y.^2)./20;
meshc(x,y,z)
 
 
五、检查求解过程（options 设置）
 
 
 
优化选项，指定为 optimoptions 的输出或 optimset 等返回的结构体。
 
 
语法：
 
options = optimoptions(@fminunc,'Algorithm','quasi-newton');
 
其中’Algorithm’设置优化工具包用什么优化算法求解问题，其中fminunc函数可以用拟牛顿法和信赖域算法。
 设置fminunc采用拟牛顿法求解：其中该函数’Algorithm’的默认值为’quasi-newton’
 
options = optimoptions(@fminunc,'Algorithm','quasi-newton');
 
 
 
quasi-newton 算法使用具有三次线搜索过程的 BFGS 拟牛顿法。这种拟牛顿法使用 BFGS公式来更新 Hessian 矩阵的逼近。您可以通过将 HessUpdate 选项设置为 ‘dfp’（并将 Algorithm 选项设置为 ‘quasi-newton’）来选择逼近逆 Hessian 矩阵的 DFP公式。您可以通过将 HessUpdate 设置为 ‘steepdesc’（并将 Algorithm 设置为 ‘quasi-newton’）来选择最陡下降法，尽管此设置通常效率不高。
 
 
设置fminunc采用信赖域算法求解：必须自己提供目标函数的梯度，不建议使用。
 
options = optimoptions(@fminunc,'Algorithm','trust-region');
 
 
 
trust-region 算法要求您在 fun 中提供梯度，并使用 optimoptions 将 SpecifyObjectiveGradient 设置为 true。此算法是一种子空间信赖域方法，基于 [2] 和 [3] 中所述的内部反射牛顿法。每次迭代都涉及使用预条件共轭梯度法 (PCG) 来近似求解大型线性方程组。
 
 
用法：
 
fun = @(x)x(1)*exp(-(x(1)^2 + x(2)^2)) + (x(1)^2 + x(2)^2)/20;
x0 = [1,2];
[x,fval,exitflag,output,grad,hessian] = fminunc(fun,x0，options)
 
其它options的参数表如下：
 
 
六、信赖域法实例（fminunc求解器trust-region设置）
 
 
 
fminunc 有两种算法：
 ‘quasi-newton’（默认值）
 ‘trust-region’
 在命令行中使用 optimoptions 设置 Algorithm 选项。
 
 
建议
 如果您的目标函数包括梯度，请使用 ‘Algorithm’ = ‘trust-region’，并将 SpecifyObjectiveGradient 选项设置为 true。
 否则，请使用 ‘Algorithm’ = ‘quasi-newton’。
 
6.1 信赖域实例1
 
考虑以下优化问题：
 
 
以下代码创建 rosenbrockwithgrad 函数，该函数包含梯度作为第二个输出。
 
function [f,g] = rosenbrockwithgrad(x)
% Calculate objective f
f = 100*(x(2) - x(1)^2)^2 + (1-x(1))^2;

if nargout > 1 % gradient required
    g = [-400*(x(2)-x(1)^2)*x(1) - 2*(1-x(1));
        200*(x(2)-x(1)^2)];
end
end
 
创建使用目标函数的梯度的选项。此外，将算法设置为 ‘trust-region’。
 
options = optimoptions('fminunc','Algorithm','trust-region','SpecifyObjectiveGradient',true);
 
将初始点设置为 [-1,2]。然后调用 fminunc。
 
x0 = [-1,2];
fun = @rosenbrockwithgrad;
x = fminunc(fun,x0,options)
 
结果
 
Local minimum found.

Optimization completed because the size of the gradient is less than
the value of the optimality tolerance.
x = 1×2

    1.0000    1.0000
 
6.2 信赖域实例2：使用问题结构体
 

 以下代码创建 rosenbrockwithgrad 函数，该函数包含梯度作为第二个输出。
 
function [f,g] = rosenbrockwithgrad(x)
% Calculate objective f
f = 100*(x(2) - x(1)^2)^2 + (1-x(1))^2;

if nargout > 1 % gradient required
    g = [-400*(x(2)-x(1)^2)*x(1)-2*(1-x(1));
        200*(x(2)-x(1)^2)];
end
end
 
将算法设置为 ‘trust-region’。
 
options = optimoptions('fminunc','Algorithm','trust-region','SpecifyObjectiveGradient',true);
 
创建包括初始点 x0 = [-1,2] 的问题结构体。有关此结构体中的必填字段，
 
problem.options = options;
problem.x0 = [-1,2];
problem.objective = @rosenbrockwithgrad;
problem.solver = 'fminunc';
 
求解
 
x = fminunc(problem)
 
结果
 
Local minimum found.

Optimization completed because the size of the gradient is less than
the value of the optimality tolerance.
x = 1×2

    1.0000    1.0000
 
原创不易，路过的各位大佬请点个赞