matlab语言与应用 06 代数方程与最优化

现代科学运算-matlab语言与应用

$\qquad \qquad$ 东北大学 http://www.icourses.cn/home/ (薛定宇)
《高等应用数学问题的MATLAB求解（第四版）》
代码可在matlab r2016b 运行。

06 代数方程与最优化问题的计算机求解

06.01 代数方程的图解法

代数方程的图解法
一元方程的图解法
一元方程的标准形式 f(x) = 0
ezplot()函数可以绘制函数曲线，其解为与实轴交点
二元方程的图解法
联立方程f(x, y) = 0 与 g(x, y) = 0
用 ezplot函数在同一坐标系下绘制出两个隐函数方程的曲线，其交点为联立方程的解

例6-1 图解法解一元方程
用图解法求解方程 $e^{-3t} \sin (4t + 2) + 4e^{-0.5t} \cos{2t} = 0.5$

ezplot('exp(-3*t)*sin(4*t+2)+4*exp(-0.5*t)*cos(2*t)-0.5', [0 5])
line([0, 5], [0, 0])
% 验证 t = 0.6738
t = 0.6738
exp(-3*t)*sin(4*t+2)+4*exp(-0.5*t)*cos(2*t)-0.5
% ans = -2.9852e-04

例6-2 图解法解二元方程
用图解法求解联立方程：
$\left \{ \begin{array}{l}x^2e^{-xy^2/2} + e^{-x/2}\sin(xy) = 0 \\ y^2\cos(y+x^2) + x^2e^{x+y} = 0 \end{array} \right.$

% 画第一个函数:
ezplot('x^2*exp(-x*y^2/2) + exp(-x/2)*sin(x*y)')
% 画第二个函数:
hold on;
ezplot('y^2*cos(y+x^2) + x^2*exp(x+y)')
hold off;

例6-3 图解法解方法
试用图解法求解二元方程
$\left \{ \begin{array}{l}x^2 + y^2 - 1 = 0 \\ 0.75x^3 -y + 0.9 = 0 \end{array} \right.$

ezplot('x^2 + y^2 - 1');
hold on;
ezplot('0.75*x^3 - y + 0.9')
hold off;
% 丢掉了2对复数根

06.02 多项式方程准解析解法

多项式方程的准解析解法
Abel-Ruffini定理说明：5阶及以上的多项式型方程没有解析解方法
特殊的高阶方程如多项式方程，可以用solve函数直接解出，必要时用vpasolve函数求准解析解
传统的数值算法得出的解不精确，这里只考虑多项式型方程的高精度求解问题

方程求解的matlab函数
求解多项式型方程的函数调用格式
最简调用方式 S = solve(eqn1, eqn2, …, eqnn)
直接得出根: [x,y, …] = solve(eqn1, eqn2, …, eqnn)
指定变量：[x, y, …] = solve(eqn1, eqn2, …, eqnn, ‘x, y, …’)

例6-4 鸡兔同笼问题
鸡兔同笼问题可以列出线性联立方程
$\left \{ \begin{array}{l}x + y = 35 \\ 2x + 4y = 94 \end{array} \right.$

syms x y;
[x0 y0] = solve(x+y == 35, 2*x+4*y==94)

二元方程问题 $\left \{ \begin{array}{l}x^2 + y^2 - 1 = 0 \\ 0.75x^3 - y + 0.9 = 0 \end{array} \right.$

[x1, y1] = vpasolve(x^2+y^2-1==0, 0.75*x^3-y+0.9==0)

例6-6 三元方程解析解
$\left \{ \begin{array}{l}x+3y^3 + 2z^2 = 1/2 \\ x^2 + 3y + z^3 = 2 \\ x^3 + 2z + 2y^2 = 2/4 \end{array} \right.$

syms x y z;
F = [x+3*y^3+2*z^2-1/2, x^2+3*y+z^3-2, x^3+2*z+2*y^2-2/4];
[x0, y0, z0] = vpasolve(F, [x, y, z])
size(x0)
% 检验
norm(subs(F, {x, y, z}, {x0, y0, z0}))

例6-7 复杂方程的准解析解
$\left \{ \begin{array}{l}\dfrac{1}{2}x^2 + x + \dfrac{3}{2} + 2\dfrac{1}{y} + \dfrac{5}{2y^2} + 3\dfrac{1}{x^3} = 0 \\ \dfrac{y}{2} + \dfrac{3}{2x} + \dfrac{1}{x^4} + 5y^4 = 0 \end{array} \right.$

syms x y;
F = [x^2/2 + x + 3/2 + 2/y + 5/(2*y^2)+3/x^3;y/2+3/(2*x)+1/x^4+5*y^4];
[x0, y0] = vpasolve(F),
size(x0),
%检验
e = norm(subs(F, {x, y}, {x0, y0}))

例6-8 含有参数的方程
$\left \{ \begin{array}{l}x^2 + ax^2 + 6b + 3y^2 = 0 \\ y = a + x + 3 \end{array} \right.$

syms x y;
syms a b;
[x1, y1] = solve(x^2+a*x^2 + 6*b + 3*y^2 == 0, y==a + x + 3, [x, y])
% [x, y] 指明要求的变量是x, y,而不是a, b

更高此带有参数的方程没有解析解

准解析解方法的优势和劣势
优势:可以求出解析解或高精度数值解，可以指定初值，可以同时求出方程的实数根和复数根
劣势:使用于多项式方程，求解其他类型方程没有优势，相比后面介绍的数值解方法速度较慢

06.03 非线性方程的数值解法

一般非线性方程数值解
求多元方程的一个实数根
s = fsolve(fun, x0)
[x, f, flag, out] = fsolve(fun, x0, opt)
[x, f, flag,out] = fsolve(fun, x0, opt, p1, p2, …)

求解代数方程组的步骤
设置变量，使等式变成如下形式
y = f(x) = 0
按如下方程描述等式
M-函数
匿名函数
inline函数，不推荐
求解方程组
验证解的正确性

例6-9 方程的数值解
$\left \{ \begin{array}{l}x^2 + y^2 - 1 = 0 \\ 0.75x^3 - y + 0.9 = 0 \end{array} \right.$
选择变量 $x_1 = x, x_2 = y$
把原始方程组变为
$\left \{ \begin{array}{l} x_1^2 + x_2^2 - 1 = 0 \\ 0.75x_1^3 - x_2 + 0.9 = 0 \end{array} \right.$
变成矩阵形式
$y = f(x) \quad y = \begin{bmatrix}x_1^2 + x_2^2 - 1 \\ 0.75x_1^3 - x_2 + 0.9 \end{bmatrix} = 0$

% m-函数
function y = myfun(x)
  y = [x(1)*x(1) + x(2)*x(2) - 1; ...
    0.75*x(1)^3 - x(2) + 0.9];
 % 匿名函数
f = @(x) [x(1)*x(1) + x(2)*x(2) - 1; ...
  0.75*x(1)^3 - x(2) + 0.9];
% inline函数
f = inline('[x(1)*x(1) + x(2)*x(2) - 1; ...
  0.75*x(1)^3 - x(2) + 0.9];', 'x');

在初值下搜索根
当初值选为 $x_0 = [1, 2]^T$

f = @(x) [x(1)*x(1) + x(2)*x(2) - 1; ...
  0.75*x(1)^3 - x(2) + 0.9];
OPT = optimset;
OPT.LargeScale = 'off';
[x, Y, c, d] = fsolve(f, [1; 2], OPT)

当使用另一个搜索初始点 $x_0 = [-1, 0]^T$

[x, Y, c, d] = fsolve(f, [-1, 0]', OPT);
x, norm(Y), kk = d.funcCount

注意：选择不同的初值可能得出不同的结果

求解方程的控制参数
选择方法和修改控制精度的函数调用格式
获得默认的常用变量
opt = optimset;
设置控制参数
直接修改: opt.TolX = 1e-10
或 set(opt, ‘TolX’, 1e-10)
其它重要参数: TolFun, MaxIter, MaxFcnEvals等

控制求解精度
重新设置相关精度的控制变量

ff = optimset;
ff.TolX = 1e-16;
ff.TolFun = 1e-30;
[x, Y, c, d] = fsolve(f, [1; 2], ff)

所期望的精度可能过于严苛，无法达到，但能尽可能精确
可以得到双精度下的最好结果

06.04 多解矩阵方程通用解法

求解多解方程的全部解
Riccati方程(第四章) $A^TX + XA - XBX + C = 0$
更多的非线性矩阵方程，例如,
广义Riccati方程 $AX + XD - XBX + C = 0$
类Riccati方程 $AX + XD - XBX^T + C = 0$
还有很多其它矩阵方程

一种简单但很使用的思路
已知:选择了初值，则fsolve函数可以求解方程，一般能得出方程的数值解
如果初值选择为复数则有可能得出复数解
可以建立一个循环，甚至是死循环
随机选择初值，调用fsolve求解
若得出的解已被记录，则放弃该解，否则记录该解
如果一段时间没有找到新的解，则终止程序
用户可以用Ctrl+C键随时终端程序
函数调用：more_sols(f, $X_0$, A, $\epsilon$, $t_{\lim}$)

   function more_sols(f,X0,varargin)
   [A,tol,tlim]=default_vals({1000,eps,30},varargin{:}); [n,m,i]=size(X0);
   if length(A)==1, a=-0.5*A; b=0.5*A; else, a=A(1); b=A(2); end
   ar=real(a); br=real(b); ai=imag(a); bi=imag(b);
   ff=optimset; ff.Display='off'; ff1=ff; ff.TolX=tol; ff.TolFun=tol; X=X0; 
   try, err=evalin('base','err'); catch, err=0; end, if i<=1; err=0; end, tic
   while (1)
      x0=ar+(br-ar)*rand(n,m);
      if abs(imag(A))>1e-5, x0=x0+(ai+(bi-ai)*rand(n,m))*1i; end
      [x,aa,key]=fsolve(f,x0,ff1);
      t=toc; if t>tlim, break; end
      if key>0, N=size(X,3);
         for j=1:N, if norm(X(:,:,j)-x)<1e-5; key=0; break; end, end
         if key>0, [x1,aa,key]=fsolve(f,x,ff);
            if norm(x-x1)<1e-5 & key>0; X(:,:,i+1)=x1; assignin('base','X',X); 
               err=max([norm(aa),err]); assignin('base','err',err); i=i+1, tic
   end, end, end, end

例6-11 求Riccati方程全部的根
求解下列Riccati方程组:
$A^TX +XA - XBX + C = 0$
其中
$A = \begin{bmatrix}-2 & 1 & -3 \\ -1 & 0 & -2 \\ 0 & -1 & -2 \end{bmatrix}, B = \begin{bmatrix}2 & 2 & -2 \\ -1 & 5 & -2 \\ -1 & 1 & 2 \end{bmatrix}, C = \begin{bmatrix}5 & -4 & 4 \\ 1 & 0 & 4 \\ 1 & -1 & 5 \end{bmatrix}$

A = [-2, 1, -3; -1, 0, -2; 0, -1, -2];
B = [2, 2, -2; -1, 5, -2; -1, 1, 2];
C = [5, -4, 4; 1, 0, 4; 1, -1, 5];
f = @(X)A'*X+X*A-X*B*X+C;
more_sols(f, zeros(3, 3, 0));
X, err

原方程有8个实根，Ctrl+C中断求解或等几分钟自动停止
求复数根

more_sols(f, X, 1000+1000i);
X, err

例6-12 变形Riccati方程求解
变形Riccati方程 $AX + XD - XBX^T +C = 0$
$A = \begin{bmatrix}2 & 1 & 9 \\ 9 & 7 & 9 \\ 6 & 5 & 3 \end{bmatrix}, B = \begin{bmatrix}0 & 3 & 6 \\ 8 & 2 & 0 \\ 8 & 2 & 8 \end{bmatrix}, C = \begin{bmatrix}7 & 0 & 3 \\ 5 & 6 & 4 \\ 1 & 4 & 4 \end{bmatrix}, D = \begin{bmatrix}3 & 9 & 5 \\ 1 & 2 & 9 \\ 3 & 3 & 0 \end{bmatrix}$
求出并检验全部根

A = [2, 1, 9; 9, 7, 9; 6, 5, 3];
B = [0, 3, 6; 8, 2, 0; 8, 2, 8];
C = [7, 0, 3; 5, 6, 4; 1, 4, 4];
D = [3, 9, 5; 1, 2, 9; 3, 3, 0];
f = @(X) A*X + X*D - X*B*X.' + C;
more_sols(f, zeros(3, 3, 0));
X, err
more_sols(f, X, 1000+1000i);
X, err

原方程有16个实根，还可以求出复数根，共40个根

例6-13 联立方程求解
求非线性方程组全部根
$\left \{ \begin{array}{l}x^2e^{-xy^2/2} + e^{-x/2} \sin(xy) = 0 \\ y^2 \cos(y+x^2) + x^2 e^{x+y} = 0 \end{array} \right.$
感兴趣区域 $-2 \pi \leq x, y \leq 2 \pi$

% 图解法
ezplot('x^2*exp(-x*y^2/2) + exp(-x/2)*sin(x*y)')
hold on;
ezplot('y^2*cos(y+x^2) + x^2*exp(x+y)')
hold off;

f = @(x) [x(1)^2*exp(-x(1)*x(2)^2/2) + ...
  exp(-x(1)/2)*sin(x(1)*x(2));
  x(2)^2*cos(x(2)+x(1)^2) + ...
  x(1)^2 * exp(x(1) + x(2))];
more_sols(f, [0; 0], [-2*pi, 2*pi]);
err

用图解法显示找出的全部根

x = X(1, 1, :);
x = x(:);
y = X(2, 1, :);
y = y(:);
plot(x, y, 'o')

高精度求解函数
将more_sols中的fsolve用vpasolve取代

  function more_vpasols(f,X0,varargin)
  [A,tlim]=default_vals({1000,60},varargin{:}); X=X0;
  if length(A)==1, a=-0.5*A; b=0.5*A; else, a=A(1); b=A(2); end
  ar=real(a); br=real(b); ai=imag(a); bi=imag(b); [i,n]=size(X0); tic
  while (1),
     x0=ar+(br-ar)*rand(1,n);
     if abs(imag(A))>1e-5, x0=x0+(ai+(bi-ai)*rand(1,n))*1i; end
     V=vpasolve(f,x0); N=size(X,1); key=1; x=sol2vec(V); 
     if length(x)>0
        t=toc; if t>tlim, break; end
        for j=1:N, if norm(X(j,:)-x)<1e-5; key=0; break; end, end
        if key>0, i=i+1;
           X=[X; x]; disp(['i=',int2str(i)]); assignin('base','X',X); tic 
  end, end, end
  function v=sol2vec(A)
  v=[]; A=struct2cell(A); for i=1:length(A), v=[v, A{i}]; end

06.05 无约束最优化问题

无约束最优化问题的数学描述
无约束最小化问题的数学描述 $\min \limits_{x} f(x)$
目标函数是一个标量函数 f(.)
向量 $x = [x_1, x_2, \cdots, x_n]^T$
决策变量，或优化变量
物理意义：求取一组x向量，使得最优化目标函数f(x)为最小
最大化问题: $\max \limits_x f(x) = \min \limits_x - f(x)$

解析解法和图解法
无约束最优化问题的必要条件：
$\dfrac{\partial f}{\partial x_1} |_{x = x^{*}} = 0, \dfrac{\partial f}{\partial x_2} |_{x - x^{*}} = 0, \cdots, \dfrac{\partial f}{\partial x_n} |_{x = x^{*}} = 0$
其中，$x^{*}$ 是最优点
方程的求解可能比解方程更麻烦，有时可能需要二阶导数运算

基于matlab的数值解法
数值最优化函数调用格式
x = fminunc(fun, x0)
x = fminsearch(fun, x0)
[x, f, flag, out] = fminsearch(fun, x0, opt, p1, p2, …)
[x, f, flag, out] = fminunc(fun, x0, opt, p1, p2, …)

目标函数的三种描述方法
M-函数(入口)
function y = fun(x)
function y = fun(x, p1, p2, …)
匿名函数
f = @(x) …
f = @(x) (x, p1, p2, …) …
inline函数(不推荐使用)
在匿名函数或inline函数中无法使用中间变量

例6-21 简单最优化问题举例
目标函数: $z = (x^2 - 2x)e^{-x^2 - y^2 -xy}$
变量替换: $x_1 = x, x_2 = y$
新目标函数: $f(x) = (x_1^2 - 2x_1)e^{-x_1^2-x_2^2-x_1x_2}$

f = @(x) (x(1)^2-2*x(1))*exp(-x(1)^2-x(2)^2-x(1)*x(2));
x0 = [0; 0];
[x, b, c, d] = fminsearch(f, x0)
% 使用fminunc()
[x, b, c, d] = fminunc(f, x0)

中间过程提取
编写一个提取、绘制中间点的函数myout,将属性Output设置成@myout

function stop = myout(x, optimValues, state)
  stop = false;
  switch state
    case 'init', hold on;
    case 'iter', plot(x(1), x(2), 'o');
      text(x(1)+0.1, x(2), int2str(optimValues.interation));
    case 'done', hold off;
  end;

搜索中间结果：
绘制三维等高线获得并叠印中间结果

f = @(x) (x(1)^2-2*x(1))*exp(-x(1)^2-x(2)^2-x(1)*x(2));
[x, y] = meshgrid(-3: .1: 3, -2: .1: 2);
z = (x.^2-2*x).*exp(-x.^2-y.^2-x.*y);
contour(x, y, z, 30);
ff = optimset;
ff.OutputFcn = @myout;
x0 = [2 1];
x = fminunc(f, x0, ff)

最优化求解函数的另一种调用方法
建立最优化问题的“结构体”模型

例6-22 用结构体模型求解

problem.solver = 'fminunc';
problem.options = optimset;
problem.objective = @(x)(x(1)^2-2*x(1))*exp(-x(1)^2-x(2)^2-x(1)*x(2));
problem.x0 = [2; 1];
[x, b, c, d] = fminunc(problem)

06.06 全局最优解

全局最优解与局部最优解
最小值存在的必要条件是 $\dfrac{d f(x)}{d x} = 0$
使用搜索方法，从初始值出发，可能找到一个这样的点，它是局部最小值
局部极小值中目标函数最小的为全局最小
整个目标函数可能存在多个局部最小值
搜索算法不一定能求出全局最小值

例6-23 基准测试函数 – Rastrigin函数
目标函数: $f(x_1, x_2) = 20 + x_1^2 + x_2^2 - 10(\cos \pi x_1 + \cos \pi x_2)$
目标函数曲面绘制：

ezsurf('20 + x1^2+x2^2 - 10*(cos(pi*x1) + cos(pi*x2))')
% view(0, 90), shading flat % 俯视图
% view(0, 0), shading flat % 侧视图

选择不同初值，得出不同结果

一个全局最优解搜索函数
[x, $f_{\min}$] = fminunc_global(fun, a, b, n, N)

function [x, f0] = fminunc_global(f, a, b, n, N, varargin)
  k0 = 0;
  f0 = Inf;
  if strcmp(class(f), 'struct')
    k0 = 1;
  end;
  for i = 1: N
    x0 = a+(b-a)*rand(n, 1);
    if k0 == 1
      f.x0 = x0;
      [x1 f1 key] = fminunc(f);
    else
      [x1 f1 key] = fminunc(f, x0, varargin{:});
    end;
    if key > 0 & f1 < f0
      x = x1;
      f0 = f1;
    end;
  end;

例6-25 改进的Rastrigin函数
将全局最优值做一个偏移
$f(x_1, x_2) = 20 + (x_1/30 -1)^2 + (x_2/20 - 1)^2 - 10[\cos(x_1/30 - 1)\pi + \cos (x_2/20 -1)\pi]$
搜索100次，97%都能找到全局最优解,耗时1分钟

f = @(x) 20 + (x(1)/30 -1)^2 + (x(2)/20-1)^2 - ...
  10*(cos(pi*(x(1)/30-1)) + cos(pi*(x(2)/20-1)));
F = [];
tic,
for i = 1:100
  [x, f0] = fminunc_global(f, -100, 100, 2, 50);
  F = [F, f0];
end;
toc

利用梯度求解最优化问题
有时，仅利用目标函数提供信息，很难得到精确的最优解。这是由于求解某些最优化问题收敛速度一般较慢，尤其是变量较多的最优化问题
可以利用梯度信息解决上述问题

例6-26 Rosenbrock函数
求Rosenbrock函数的无约束最优化问题
人造函数:$f(x_1, x_2) = 100(x_2 - x_1^2)^2 + (1 - x_1)^2$

% 绘制三维等高线图(香蕉函数)
[x, y] = meshgrid(0.5: 0.01: 1.5);
z = 100*(y - x.^2).^2 +(1-x).^2;
contour3(x, y, z, 100),
zlim([0, 310])
% 不用梯度信息
f = @(x)100*(x(2)-x(1)^2)^2 +(1-x(1))^2;
ff = optimset;
ff.TolX = 1e-10;
ff.TolFun = 1e-20;
x = fminunc(f, [0; 0], ff)

采用梯度信息求解
求梯度矩阵

syms x1 x2;
f = 100*(x2-x1^2)^2 + (1 - x1)^2;
J = jacobian(f, [x1, x2])

自包含梯度的目标函数

function [y, Gy] = c6fun3(x)
  y = 100*(x(2)-x(1)^2)^2 + (1-x(1))^2;
  Gy = [-400*(x(2)-x(1)^2)*x(1)-2+2*x(1);
    200*x(2)-200*x(1)^2];

求解

ff.GradObj = 'on';
x = fminunc(@c6fun3, [0, 0], ff)

06.07 可行解区域

有约束最优化问题的计算机求解

约束条件与可行解区域
有约束非线性最优化问题的一般描述
$\min \limits_{x \; s.t. \; G(x) \leq 0} f(x)$
决策变量： $x = [x_1, x_2, \cdots, x_n]^T$
约束条件：$G(x) \leq 0$
满足约束条件的所有x称为可行解区域
在满足约束条件($G(x \leq 0$)的前提下最优

例6-28 图解法求最优化
图解法 $\max \limits_{x \; s.t. \; \left \{ \begin{array}{l} 9 \geq x_1^2 + x_2^2 \\ x_1 + x_2 \leq 1 \end{array} \right.} (-x_1^2 -x_2)$
目标函数描述：

[x1, x2] = meshgrid(-3: .1: 3);
z = -x1.^2 - x2;

可行解区域描述

i = find(x1.^2 + x2.^2 > 9);
z(i) = NaN;
i = find(x1 + x2 > 1);
z(i) = NaN;
surf(x1, x2, z);
shading interp;

带有变量边界约束的最优化问题求解
带有变量边界约束的最优化问题的数学描述
$\min \limits_{x \; s.t. \; x_m \leq x \leq x_M} f(x)$
其中，符号s.t.表示subject to
fminsearchbnd()函数在本书程序包中给出或从Mathworks的File Exchange下载

function [x,fval,exitflag,output]=fminsearchbnd3(fun,x0,LB,UB,options,varargin)
xsize = size(x0);
x0 = x0(:);
n=length(x0);

if (nargin<3) || isempty(LB)
  LB = repmat(-inf,n,1);
else
  LB = LB(:);
end
if (nargin<4) || isempty(UB)
  UB = repmat(inf,n,1);
else
  UB = UB(:);
end

if (n~=length(LB)) || (n~=length(UB))
  error 'x0 is incompatible in size with either LB or UB.'
end

% set default options if necessary
if (nargin<5) || isempty(options)
  options = optimset('fminsearch');
end

% stuff into a struct to pass around
params.args = varargin;
params.LB = LB;
params.UB = UB;
params.fun = fun;
params.n = n;
params.OutputFcn = [];

params.BoundClass = zeros(n,1);
for i=1:n
  k = isfinite(LB(i)) + 2*isfinite(UB(i));
  params.BoundClass(i) = k;
  if (k==3) && (LB(i)==UB(i))
    params.BoundClass(i) = 4;
  end
end

x0u = x0;
k=1;
for i = 1:n
  switch params.BoundClass(i)
    case 1
      % lower bound only
      if x0(i)<=LB(i)
        % infeasible starting value. Use bound.
        x0u(k) = 0;
      else
        x0u(k) = sqrt(x0(i) - LB(i));
      end

      % increment k
      k=k+1;
    case 2
      % upper bound only
      if x0(i)>=UB(i)
        % infeasible starting value. use bound.
        x0u(k) = 0;
      else
        x0u(k) = sqrt(UB(i) - x0(i));
      end

      % increment k
      k=k+1;
    case 3
      % lower and upper bounds
      if x0(i)<=LB(i)
        % infeasible starting value
        x0u(k) = -pi/2;
      elseif x0(i)>=UB(i)
        % infeasible starting value
        x0u(k) = pi/2;
      else
        x0u(k) = 2*(x0(i) - LB(i))/(UB(i)-LB(i)) - 1;
        % shift by 2*pi to avoid problems at zero in fminsearch
        % otherwise, the initial simplex is vanishingly small
        x0u(k) = 2*pi+asin(max(-1,min(1,x0u(k))));
      end

      % increment k
      k=k+1;
    case 0
      % unconstrained variable. x0u(i) is set.
      x0u(k) = x0(i);

      % increment k
      k=k+1;
    case 4
      % fixed variable. drop it before fminsearch sees it.
      % k is not incremented for this variable.
  end

end
% if any of the unknowns were fixed, then we need to shorten
% x0u now.
if k<=n
  x0u(k:n) = [];
end

% were all the variables fixed?
if isempty(x0u)
  % All variables were fixed. quit immediately, setting the
  % appropriate parameters, then return.

  % undo the variable transformations into the original space
  x = xtransform(x0u,params);

  % final reshape
  x = reshape(x,xsize);

  % stuff fval with the final value
  fval = feval(params.fun,x,params.args{:});

  % fminsearchbnd was not called
  exitflag = 0;

  output.iterations = 0;
  output.funcount = 1;
  output.algorithm = 'fminsearch';
  output.message = 'All variables were held fixed by the applied bounds';

  % return with no call at all to fminsearch
  return
end

% Check for an outputfcn. If there is any, then substitute my
% own wrapper function.
if ~isempty(options.OutputFcn)
  params.OutputFcn = options.OutputFcn;
  options.OutputFcn = @outfun_wrapper;
end

% now we can call fminsearch, but with our own
% intra-objective function.
[xu,fval,exitflag,output] = fminsearch(@intrafun,x0u,options,params);

% undo the variable transformations into the original space
x = xtransform(xu,params);

% final reshape
x = reshape(x,xsize);

% Use a nested function as the OutputFcn wrapper
  function stop = outfun_wrapper(x,varargin);
    % we need to transform x first
    xtrans = xtransform(x,params);

    % then call the user supplied OutputFcn
    stop = params.OutputFcn(xtrans,varargin{1:(end-1)});

  end

end % mainline end

% ======================================
% ========= begin subfunctions =========
% ======================================
function fval = intrafun(x,params)
% transform variables, then call original function

% transform
xtrans = xtransform(x,params);

% and call fun
fval = feval(params.fun,xtrans,params.args{:});

end % sub function intrafun end

% ======================================
function xtrans = xtransform(x,params)
% converts unconstrained variables into their original domains

xtrans = zeros(1,params.n);
% k allows some variables to be fixed, thus dropped from the
% optimization.
k=1;
for i = 1:params.n
  switch params.BoundClass(i)
    case 1
      % lower bound only
      xtrans(i) = params.LB(i) + x(k).^2;

      k=k+1;
    case 2
      % upper bound only
      xtrans(i) = params.UB(i) - x(k).^2;

      k=k+1;
    case 3
      % lower and upper bounds
      xtrans(i) = (sin(x(k))+1)/2;
      xtrans(i) = xtrans(i)*(params.UB(i) - params.LB(i)) + params.LB(i);
      % just in case of any floating point problems
      xtrans(i) = max(params.LB(i),min(params.UB(i),xtrans(i)));

      k=k+1;
    case 4
      % fixed variable, bounds are equal, set it at either bound
      xtrans(i) = params.LB(i);
    case 0
      % unconstrained variable.
      xtrans(i) = x(k);

      k=k+1;
  end
end

end % sub function xtransform end

matlab求解方法
变量边界约束的最优化问题求解的语句调用格式
x = fminsearchbnd(fun, x0, xm, xM)
[x, f, flag, out] = fminsearchbnd(fun, x0, xm, xM, opt, p1, p2, …)

例6-29 边界约束问题求解
求解 Resenbrock问题
$f(x_1, x_2) = 100(x_2 - x_1^2)^2 + (1 - x_1)^2$
其中，$x_1 \in (2, 4)$ 和 $x_2 \in (3, 6)$

f = @(x)[100*(x(2)-x(1)^2)^2 + (1-x(1))^2];
xm = [2, 3];
xM = [4, 6];
x = fminsearchbnd3(f, [0, 0], xm, xM)

06.08 线性规划与二次型规划

线性规划问题的计算机求解
线性规划(LP)问题的一般数学描述为:
$\min \limits _{x \; s.t. \; \left \{ \begin{array}{l} Ax \leq B \\ A_{eq} x = B_{eq} \\ x_m \leq x \leq x_M \end{array} \right.} f^T x$
目标函数和约束条件都是线性的
注意，约束的标准形式是 $\leq$
线性规划是凸问题，与初值无关
求解线性规划问题的函数调用：
x = linrog(f, A, B)
x = linprog(f, A, B, $A_{eq}$, $B_{eq}$)
x = linprog(f, A, B, $A_{eq}, B_{eq}, x_m, x_M, x_0$)
[x, $f_{opt}$, flag, c] = linprog(f, A, B, $A_{eq}, B_{eq}, x_m, x_M, x_0$, OPT)
x = linprog(problem)
[x, $f_{opt}$, flag, c] = linprog(problem)

例6-29 线性规划问题求解
试求解下面线性规划问题
$\min \limits_{x \; s.t. \; \left \{ \begin{array}{l}2x_2 + x_3 + 4x_4 + 2x_5 \leq 54 \\ 3x_1 + 4x_2 + 5x_3 -x_4 - x_5 \leq 64 \\ x_1, x_2 \geq 0, x_3 \geq 3.32, x_4 \geq 0.678, x_5 \geq 2.57 \end{array} \right.} (-2x_1 -x_2 - 4x_3 - 3x_4 -x_5)$
约束条件的矩阵表示: $f = [-2, -1, -4, -3, -1]^T$
$A = \begin{bmatrix}0 & 2 & 1 & 4 & 2 \\ 3 & 4 & 5 & -1 & -1 \end{bmatrix}, B = \begin{bmatrix}54 \\ 62 \end{bmatrix}$

f = -[2, 1, 4, 3, 1]';
Ae = [];
Be = [];
A = [0, 2, 1, 4, 2; 3,4,5, -1, -1];
B = [54; 62];
xm = [0, 0, 3.32, 0.678, 2.57];
[x, f_opt, key, c] = linprog(f, A, B, Ae, Be, xm, [], [])

例 6-30 结构体求解

clear P;
P.f = f;
P.Aineq = A;
P.Bineq = B;
P.lb = xm;
P.solver = 'linprog';
P.options = optimset;
[x, f_opt, key, c] = linprog(P)

例6-31 线性规划问题求解
求解线性规划问题
$\max \limits_{x \; s.t. \; \left \{ \begin{array}{l} x_1/4 - 60 x_2 - x_3/50 + 9x_4 \leq 0 \\ -x_1/2 + 90 x_2 + x_5/50 - 3x_4 \geq 0 \\ x_3 \leq 1, x_1 \geq -5, x_2 \geq -5, x_3 \geq -5, x_4 \geq -5 \end{array} \right.} [\dfrac{3}{4} x_1 - 150 x_2 + \dfrac{1}{50} x_3 - 6x_4]$
先将原问题转换为标准问题
$f = \begin{bmatrix} -\dfrac{3}{4} & 150 & -\dfrac{1}{50} & 6 \end{bmatrix}$
$A = \begin{bmatrix} 0.25 & -60 & -\dfrac{1}{50} & 9 \\ 0.5 & -90 & -\dfrac{1}{50} & 3 \end{bmatrix}, B = \begin{bmatrix} 0 \\ 0 \end{bmatrix}, x_M = \begin{bmatrix} \infty \\ \infty \\ 1 \\ \infty \end{bmatrix}$

f = [-3/4, 150, -1/50, 6];
Aeq = [];
Beq = [];
A = [1/4, -60, -1/50, 9; 1/2, -90, -1/50, 3];
B = [0, 0];
xm = [-5; -5; -5; -5];
xM = [Inf; Inf; 1; Inf];
[x, f_opt, key, c] = linprog(f, A, B, Aeq, Beq, xm, xM, [0; 0; 0; 0])

% 结构体求解
clear P;
P.f = f;
P.Aineq = A;
P.Bineq = B;
P.solver = 'linprog';
P.lb = xm;
P.ub = xM;
P.options = optimset;
[x, f_opt, key, c] = linprog(P)

例6-32 双下标问题
求解双下标线性规划问题
$\min \quad 2800(x_{11} + x_{21} + x_{31} + x_{41}) + 4500(x_{12} + x_{22} + x_{32}) + 6000(x_{13} + x_{23}) + 7300x_{14}$
$x \; s.t. \; \left \{ \begin{array}{l}x_{11} + x_{12} + x_{13} + x_{14} \geq 15 \\ x_{12} + x_{13} + x_{14} + x_{21} + x_{22} + x_{23} \geq 10 \\ x_{13} + x_{14} + x_{22} + x_{23} + x_{31} + x_{32} \geq 20 \\ x_{14} + x_{23} + x_{32} + x_{41} \geq 12 \\ x_{ij} \geq 0, (i = 1, 2, 3, 4, j = 1, 2, 3, 4) \end{array} \right.$
变成标准型，定义为单下标变量
$x_1 = x_{11}, x_2 = x_{12}, x_3 = x_{13}, x_4 = x_{14}, x_5 = x_{21},$
$x_6 = x_{22}, x_7 = x_{23}, x_8 = x_{31}, x_9 = x_{32}, x_10 = x_{41}$
原问题改为：
$\min \quad 2800(x_1 + x_5 + x_8 + x_{10}) + 4500(x_2 + x_6 + x_9) + 6000(x_3 + x_7) + 7300x_4$
$x \; s.t. \; \left \{ \begin{array}{l}-（x_1 + x_2 + x_3 + x_4) \leq 15 \\ -(x_2 + x_3 + x_4 + x_5 + x_6 + x_7) \leq 10 \\ -(x_3 + x_4 + x_6 + x_7 + x_8 + x_9) \leq 20 \\ -(x_4 + x_7 + x_9 + x_{10}) \leq 12 \\ x_i \geq 0, (i = 1, 2, 3, \cdots, 10) \end{array} \right.$

f = 2800*[1 0 0 0 1 0 0 1 0 1] + ...
     4500*[0 1 0 0 0 1 0 0 1 0] + ...
     6000*[0 0 1 0 0 0 1 0 0 0] + ...
     7300*[0 0 0 1 0 0 0 0 0 0];
A = -[1 1 1 1 0 0 0 0 0 0; 0 1 1 1 1 1 1 0 0 0;
         0 0 1 1 0 1 1 1 1 0; 0 0 0 1 0 0 1 0 1 1];
B = -[15; 10; 20; 12];
xm = [0 0 0 0 0 0 0 0 0 0];
Aeq = [];
Beq = [];
x = linprog(f, A, B, Aeq, Beq, xm)

得出的结果再反代回双下标自变量

二次型规划求解
一般二次型绘画问题的数学表示
$\min \limits_{x \; s.t. \; \left \{ \begin{array}{l}Ax \leq B \\ A_{eq} x = B_{eq} \\ x_m \leq x \leq x_M \end{array} \right.} \Big(\dfrac{1}{2} x^T Hx + f^T x \Big)$
二次型规划问题也是凸问题，其解与初值选择无关
x = quadprog(H, f, A, B)
x = quadprog(H, f, A, B, $A_{eq}, B_{eq}$)
x = quadprog(H, f, A, B, $A_{eq}, B_{eq}, x_m, x_M$)
[x, $f_{opt}$, flag, c] = quadprog(H, f, A, B, $A_{eq}, B_{eq}, x_m, x_M, x_0$, OPT)
x = quadprog(progblem)
[x, $f_{opt}$, flag, c] = quadprog(progblem)
构造H矩阵注意标准型系数1/2

06.09 非线性规划

一般非线性规划问题
一般非线性规划问题
$\min \limits_{x \; s.t. \; G(x) \leq 0} f(x)$
其中：$x = [x_1, x_2, \cdots, x_n]^T$
物理解释：在给出的约束条件下，找出向量x，使目标函数达到最小值
决策变量
可行解
更具体描述
$\min \limits_{x \; s.t. \; \left \{ \begin{array}{l}Ax \leq B \\ A_{eq}x = B_{eq} \\ x_m \leq x \leq x_M \\ C(x) \leq 0 \\ C_{eq}(x) = 0 \end{array} \right.} f(x)$
求解出非线性规划问题：
[x, $f_{opt}$, flag, c] = fmincon(fun, $x_0, A, B, A_{eq}, B_{eq}, x_m, x_M, CFun, OPT, p_1, p_2, \cdots$)
CFun返回两个输出变量，不等式约束C和等式约束$C_{eq}$,不能用匿名函数来描述，需要使用
.m文件

例6-34 非线性规划
试求解下面非线性规划问题
$\min \quad 1000 - x_1^2 - 2x_2^2 - x_3^2 - x_1x_2 - x_1x_3$
$x \; s.t. \; \left \{ \begin{array}{l}x_1^2 + x_2^2 +x_3^2 - 25 = 0 \\ 8x_1 + 14x_2 + 7x_3 - 56 = 0 \\ x_1, x_2, x_3 \geq 0 \end{array} \right.$
为目标函数和约束函数编写M-函数
目标函数

% 非线性约束函数(返回两个变量)
function [c, ceq] = opt_con1(x)
  ceq = [x(1)*x(1) + x(2)*x(2)+x(3)*x(3)-25;
    8*x(1) + 14*x(2) + 7*x(3) - 56];
  c = [];

% 目标函数
f = @(x) 1000-x(1)*x(1) - 2*x(2)*x(2) - x(3)*x(3) - x(1)*x(2) -x(1)*x(3);
% 问题求解
ff = optimset;
ff.LargeScale = 'off';
ff.TolFun = 1e-30;
ff.TolX = 1e-15;
ff.TolCon = 1e-20;
x0 = [1;1;1];
xm = [0; 0; 0];
xM = [];
A = [];
B = [];
Aeq = [];
Beq = [];
[x, f_opt, c, d] = fmincon(f, x0, A, B, Aeq, Beq, xm, xM, @opt_con1, ff)

% 使用结构体描述
clear P;
P.objective = f;
P.nonlcon = @opt_con1;
P.x0 = x0;
P.lb = xm;
P.options = ff;
P.solver = 'fmincon';
[x, f_opt, c, d] = fmincon(P)

分离掉线性约束条件

function [c, ceq] = opt_con2(x)
  c = [];
  ceq = x(1)*x(1) + x(2)*x(2) + x(3)*x(3) - 25;

求出结果：

x0 = [1;1;1];
Aeq = [8, 14, 7];
Beq = 56;
[x, f_opt, c, d] = fmincon(f, x0, A, B, Aeq, Beq, xm, xM, @opt_con2, ff)

例6-38 复杂优化问题
求解最优化问题
$\min \quad k$
$q, w, k \; s.t. \; \left \{ \begin{array}{l}q_3 + 9.625q_1 w + 16 q_2 w + 16 w^2 + 12 - 4q_1 - q_2 - 78 w = 0 \\ 16q_1 w + 44 - 19 q_1 - 8 q_2 - q_3 - 24 w = 0 \\ 2.25 - 0.25k \leq q_1 \leq 2.25 + 0.25k \\ 1.5 - 0.5k \leq q_2 \leq 1.5 + 0.5k \\ 1.5 - 1.5k \leq q_3 \leq 1.5 + 1.5k \end{array} \right.$
选择决策变量
$x_1 = q_1, x_2 = q_2, x_3 = q_3, x_4 = w, x_5 = k$
将原问题修改为标准型
$\min \quad x_5$
$x \; s.t. \; \left \{ \begin{array}{l}x_3 + 9.625x_1 x_4 + 16 x_2 x_4 + 16 x_4^2 + 12 - 4x_1 - x_2 - 78 x_4 = 0 \\ 16x_1 x_4 + 44 - 19 x_1 - 8 x_2 - x_3 - 24 x_4 = 0 \\ -0.25x_5 - x_1 \leq -2.25 \\ x_1 - 0.25x_5 \leq 2.25 \\ -0.5 x_5 - x_2 \leq -1.5 \\ x_2 - 0.5x_5 \leq 1.5 \\ -1.5x_5 - x_3 \leq -1.5 \\ x_3 - 1.5x_5 \leq 1.5 \end{array} \right.$

例6-39 全局最优解
例6-38 容易陷入局部最优解

function [c,ceq] = c6exnls(x)
  c = []; 
  ceq = [x(3)+9.625*x(1)*x(4)+16*x(2)*x(4)+16*x(4)^2+12-4*x(1)-x(2)-78*x(4); 
    16*x(1)*x(4)+44-19*x(1)-8*x(2)-x(3)-24*x(4)];

% 全局最优化求解函数
function [x,f0] = fmincon_global(f, a, b, n, N, varargin)
  x0 = rand(n,1);
  k0 = 0;
  if strcmp(class(f), 'struct')
    k0 = 1;
  end %处理结构体
  if k0 == 1
    f.x0 = x0;
    [x f0] = fmincon(f); %如果是结构体描述的问题，直接求解
  else
    [x f0] = fmincon(f, x0, varargin{:});
  end %如果不是结构体描述的，直接求解
  for i = 1:N
    x0 = a+(b-a)*rand(n,1); %用循环结构尝试不同的随机搜索初值
    if k0 == 1
      f.x0 = x0;
      [x1 f1 key] = fmincon(f); %结构体问题求解
    else
      [x1 f1 key] = fmincon(f, x0, varargin{:});
    end %非结构体问题求解
    if key>0 & f1<f0
      x = x1;
      f0 = f1;
    end %如果找到的解优于现有的最好解，存储该解
  end

clear P;
P.objective = @(x)x(5);
P.nonlcon = @c6exnls;
P.solver = 'fmincon';
P.Aineq = [-1, 0, 0, 0, -0.25; 1, 0, 0, 0, -0.25; 0, -1, 0, 0, -0.5;
    0, 1, 0, 0, -0.5; 0, 0, -1, 0, -1.5; 0, 0, 1, 0, -1.5];
P.Bineq = [-2.25; 2.25; -1.5; 1.5; -1.5; 1.5];
P.options = optimset;
P.x0 = rand(5, 1);
[x, f0] = fmincon_global(P, 0, 5, 5, 50)
% 运行100次，成功率100%，耗时4分钟
tic, X=[];
for i = 1:100
  [x, f0] = fmincon_global(P, 0, 5, 5, 50);
  X = [X; x'];
end;
toc

06.10 整数规划的穷举方法

混合整数规划问题的计算机求解
非线性规划包括了单目标规划的绝大多数最优化类型，在实际应用中还可能有其它要求，比如某决策变量为人数，则需要为整数，这样就引入了整数规划问题

整数规划问题的穷举方法
所谓穷举方法，就是将决策变量取值所有的可能都衡量一番，从满足条件的决策变量可能中找出目标函数值最小的组合，这样的解即为原始问题的全局最优解。
如果已知自变量所在区间，理论上可以考虑用穷举方法举出区间内所有的变量组合
不能用于NP问题，穷举方法只适合于极有限的小规模问题
混合整数规划问题不能采用穷举方法

例6-40 线性整数规划
求解下列线性整数规划问题，其中，各个变量全为整数
$\min \limits_{x \; s.t. \; \left \{ \begin{array}{l}2x_2 + x_3 + 4x_4 + 2x_3 \leq 54 \\ 3x_1 + 4x_2 + 5x_3 -x_4 -x_5 \leq 62 \\ x_1, x_2 \geq 0, x_3 \geq 3.32, x_4 \geq 0.678, x_5 \geq 2.57 \end{array} \right.} (-2x_1 -x_2 -4x_3 -3x_4 -x_5)$
穷举法思路
认为假定取值范围[0.25]
用 ndgrid 函数可以生成所有可能
对可行解排序，得出最优解和次优解

穷举求解

N = 25;
[x1, x2, x3, x4, x5] = ndgrid(0:N, 0:N, 4:N, 1:N, 3:N);
i = find((2*x2+x3+4*x4+2*x5 <= 54) & ...
    (3*x1 + 4*x2 + 5*x3 - x4 - x5 <= 62));
x1 = x1(i);
x2 = x2(i);
x3 = x3(i);
x4 = x4(i);
x5 = x5(i);
f = -2*x1 - x2 -4*x3 - 3*x4 -x5;
[fmin, ii] = sort(f);
in = ii(1);
x = [x1(in), x2(in), x3(in), x4(in), x5(in)]
% size(x1)

% 最优解与次最优解
L = 10;
f1 = f(1:L)';
in = ii(1:L);
x = [x1(in), x2(in), x3(in), x4(in), x5(in), f1']

穷举法的优劣
限定的区间[0, 15],不能随意拓展到此区间外，如想将25变为30，在一般的计算机配置下都实现不了，因为需要内除过大，5个变量的存储量为
$31^5 \times 5 \times 8 / 2^{20} = 1092.1MB$
除了得出最优解外，事实上还可以得出若干组合，为次最优解

例6-41 穷举方法求解
用穷举法求解非线性整数规划问题
$\min \quad x_1^2 + x_2^2 + 2x_3^2 + x_4^2 - 5x_1 - 5x_2 - 21x_3 + 7x_4$
$x \; s.t. \; \left \{ \begin{array}{l}-x_1^2 -x_2^2 -x_3^2 -x_4^2 -x_1 +x_2 -x_3 + x_4 + 8 \geq 0 \\ -x_1^2 -2x_2^2 -x_3^2 -2x_4^2 +x_1 +x_4 + 10 \geq 0 \\ -2x_1^2 -x_2^2 -x_3^2 -2x_4^2 +x_2 +x_4 + 5 \geq 0 \end{array} \right.$

N = 30;
[x1 x2 x3 x4] = ndgrid(-N:N);
ii = find((-x1.^2 -x2.^2 -x3.^2 -x4.^2 -x1 +x2 -x3 +x4 + 8 >= 0) & ...
    (-x1.^2 -2*x2.^2 -x3.^2 -2*x4.^2 + x1 + x4 + 10 >= 0) & ...
    (-2*x1.^2 -x2.^2 -x3.^2 -2*x4.^2 +x2 + x4 + 5 >= 0));
x1 = x1(ii);
x2 = x2(ii);
x3 = x3(ii);
x4 = x4(ii);
ff = x1.^2 + x2.^2 + 2*x3.^2 +x4.^2 -5*x1 - 5*x2 - 21*x3 +7*x4;
[fm, ii] = sort(ff);
k = ii(1:5);
X = [x1(k), x2(k), x3(k), x4(k)],
fm(1:5)

例6-42 离散规划问题求解
离散最优化问题
$\min \quad 2x_1^2 + x_2^2 - 16x_1 - 10x_2$
$x \; s.t. \; \left \{ \begin{array}{l}x_1^2 -6x_1 + x_2 - 11 \leq 0 \\ -x_1x_2 + 3x_2 + e^{x_1-3} - 1\leq 0 \end{array} \right.$
$x_1 是0.25的整数倍, x_2是0.1的整数倍$

[x1 x2] = meshgrid(-20:0.25:20, 3:0.1:20);
ii = find(x1.^2 - 6*x1 + x2 - 11 <= 0 & ...
    -x1.*x2 + 3*x2 + exp(x1 -3) - 1 <= 0);
x1 = x1(ii);
x2 = x2(ii);
ff = 2*x1.^2 + x2.^2 - 16*x1 - 10*x2;
[fm, ij] = sort(ff);
k = ij(1:5);
X = [x1(k) x2(k)],
fm(1:5)

06.11 混合整数规划问题

整数线性规划问题求解
整数线性规划的一般数学描述
$\min \limits_{x \; s.t. \; \left \{ \begin{array}{l}Ax \leq B \\ A_{eq} x = B_{eq} \\ x_m \leq x \leq x_M \\ \hat x \; integers \end{array} \right.} f^T x$
$其中，\hat x 为变量x的子集$
整数规划与混合整数规划

混合整数线性规划的求解函数
混合线性整数规划的求解函数
[x, f_m, key, c] = intlinprog(problem)
[x, f_m, key, c] = intlinprog($f, intcon, A, b, A_{eq}, b_{eq}, x_m, x_M$, pars)
结构体设置
problem.solver = ‘intlinprog’;
problem.options = optimoptions(‘intlinprog’)
结果精调
x(intcon) = round(x(intcon))

例6-42 混合整数线性规划
混合整数线性规划, 1, 4, 5变量要求为整数
$\min \quad (-2x_1 -x_2 -4x_3 -3x_4 -x_5)$
$x \; s.t. \; \left \{ \begin{array}{l}2x_2 + x_3 + 4x_4 + 2x_5 \leq 54 \\ 3x_1 + 4x_2 + 5x_3 -x_4 -x_5 \leq 62 \\ x_1, x_2 \geq 0, x_3 \geq 3.32, x_4 \geq 0.678, x_5 \geq 2.57 \end{array} \right.$

clear P;
P.solver = 'intlinprog';
P.options = optimoptions('intlinprog');
P.lb = [0; 0; 3.32; 0.678; 2.57];
P.f = [-2, -1, -4, -3, -1];
I = [1, 4, 5];
P.Aineq = [0, 2, 1, 4, 2; 3, 4, 5, -1, -1];
P.Bineq = [54; 62];
P.intcon = [1, 4, 5];
[x, f, a, b] = intlinprog(P);
x(I) = round(x(I))

一般非线性整数规划问题与求解
分枝定界算法是一种最常用的处理非线性整数规划或混合规划问题的算法
免费工具箱由 Mr.Keort Kuipers 提供
数学描述
$\min \quad f(x)$
$x \; s.t. \; \left \{ \begin{array}{l} Ax \leq B \\ A_{eq} x = B_{eq} \\ x_m \leq x \leq x_M \\ C(x) \leq 0 \\ C_{eq}(x) = 0 \\ \hat x \; integers \end{array} \right.$
求解语句：
[err, f, x] = bnb20_new(fun, $x_0， intcon, x_m, x_M, A, B, A_{eq}, B_{eq}$, CFun)

例6-46 非线性整数规划
整数规划
$\min \quad x_1^3 + x_2^2 - 4x_1 + 4 + x_3^4$
$x \; s.t. \; \left \{ \begin{array}{l}x_1 - 2x_2 + 12 + x_3 \geq 0 \\ -x_1^2 + 3x_2 - 8 - x_3 \geq 0 \\ x_1 \geq 0, x_2 \geq 0, x_3 \geq 0 \end{array} \right.$
非线性约束条件

function [c, ce] = c6exinl(x)
  ce = [];
  c = [-x(1) + 2*x(2) - 12 - x(3);
      x(1)^2-3*x(2)+8+x(3)];

clear P;
in = [1:3];
P.objective = @(x)x(1)^3+x(2)^2-4*x(1) + 4 + x(3)^4;
P.intcon = in;
P.nonlcon = @c6exinl;
P.lb = [0; 0; 0];
P.ub = 100*[1;1;1];
P.x0 = P.ub;
[err, fm x] = BNB20_new(P);
x(in) = round(x(in))

例6-47 离散规划问题求解
离散最优化问题
$\min \quad 2x_1^2 + x_2^2 - 16x_1 - 10x_2$
$x \; s.t. \; \left \{ \begin{array}{l}x_1^2 - 6x_1 +x_2 - 11 \leq 0 \\ -x_1 x_2 + 3x_2 + e^{x_1 -3 } - 1\leq 0 \end{array} \right.$
$x_1 是0.25 的整数倍， x_2 是 0.1 的整数倍$
引入 $y_1 = 4x_1, y_2 = 10x_2,$整数规划
$\min \quad 2y_1^2/16 + y_2^2/100 - 4y_1 - y_2$
$y \; s.t. \; \left \{ \begin{array}{l} y_1^2/1y - 6y_1/4 + y_2/10 - 11 \leq 0 \\ -y_1y_2/40 + 3y_2/10 + e^{y_1/4-3} -1 \leq 0 \\ y_2 \geq 30 \end{array} \right.$
约束条件

function [c, ceq] = c6mdisp(y)
    ceq = [];
    c = [y(1)^2/10-6*y(1)/4+y(2)/10-11;
        -y(1)*y(2)/40+3*y(2)/10+exp(y(1)/4-3)-1];

clear P;
P.objective = @(y)2*y(1)^2/16+y(2)^2/100 - ...
    4*y(1) - y(2);
P.nonlcon = @c6mdisp;
P.lb = [-200; 30];
P.ub = [200; 200];
P.intcon = [1; 2];
P.x0 = [12; 30];
[errmsg, ym, y] = BNB20_new(P);
x = [y(1)/4, y(2)/10]

0-1规划问题求解
没有现成的函数可以直接求解0-1规划
可以利用前面介绍的整数规划问题求解函数
线性规划问题：intlinprog
非线性规划问题: BNB20_new
使用技巧:
将下界定义一个全为0的向量
将上界定义一个全为1的向量

例6-50 混合非线性0-1规划
混合0-1规划问题
$\min \quad 5y_1 + 6y_2 + 8y_3 + 10x_1 - 7x_3 - 18 \ln (x_2+1) - 19.2 \ln(x_1 -x_2 + 1) + 10$
$x \; s.t. \; \left \{ \begin{array}{l} 0.8 \ln(x_2 + 1) + 0.96 \ln (x_1 - x_2 + 1) - 0.8 x_3 \leq 0 \\ \ln (x_2 + 1) + 1.2 \ln (x_1 - x_2 + 1) - x_3 - 2y_3 \geq -2 \\ x_2 - x_1 \leq 0 \\ x_2 - 2y_1 \leq 0 \\ x_1 - x_2 - 2y_2 \leq 0 \\ y_1 + y_2 \leq 1 \\ 0 \leq x \leq [2, 2, 1]^T, y \in \lbrace 0, 1 \rbrace \end{array} \right.$
变量替换，保留前三个x
$令 x_4 = y_1, x_5 = y_2, x_6 = y_3$
$\min \quad 5x_4 + 6x_5 + 8x_6 + 10 x_1 - 7x_3 - 18 \ln (x_2 + 1) - 19.2 \ln (x_1 - x_2 + 1) + 10$
$x \; s.t. \; \left \{ \begin{array}{l} 0.8 \ln(x_2 + 1) + 0.96 \ln (x_1 - x_2 + 1) - 0.8 x_3 \leq 0 \\ \ln (x_2 + 1) + 1.2 \ln (x_1 - x_2 + 1) - x_3 - 2x_6 \geq -2 \\ x_2 - x_1 \leq 0 \\ x_2 - 2x_4 \leq 0 \\ x_1 - x_2 - 2x_5 \leq 0 \\ x_4 + x_5 \leq 1 \\ 0 \leq x \leq [2, 2, 1, 1, 1, 1]^T \end{array} \right.$

function [c, ceq] = c6mmibp(x)
    ceq = [];
    c = [-0.8*log(x(2)+1)-0.96*log(x(1)-x(2)+1)+0.8*x(3);
        -log(x(2)+1)-1.2*log(x(1)-x(2)+1)+x(3)+2*x(6)-2];

clear P;
P.intcon = 4:6;
P.x0 = [0 0 0 0 0 0]';
P.objective = @(x)5*x(4)+6*x(5)+8*x(6)+10*x(1) ...
    -7*x(3)-18*log(x(2)+1)-19.2*log(x(1)-x(2)+1)+10;
P.ub = [2 2 1 1 1 1 ]';
P.lb = [0 0 0 0 0 0]';
P.Aineq = [-1 1 0 0 0 0; 0 1 0 -2 0 0; 1 -1 0 0 -2 0; 0 0 0 1 1 0];
P.nonlcon = @c6mmibp;
P.Bineq = [0; 0; 0; 1];
[errmsg, fm, x] = BNB20_new(P)

06.12 动态规划与最优路径问题

多阶段目标–目标函数随着实际完成情况动态地改变
图的矩阵表示方法
有向图的路径寻优
无向图的路径最优搜索

例6-57 最短路径问题
有向图:节点、边、权值

倒叙求解
问题：
繁琐、易出错、难以求解大规模问题

图的矩阵表示方法
图是由节点和边构成的。边是连接两个节点的直接路径。如果边是又向的，则图称为有向图，否则称为无向图
在计算机中，图用矩阵表示，即关联矩阵
matlab语言支持关联矩阵的稀疏矩阵表示方法

例5-58 有向图的矩阵描述

关联矩阵表示
构造关联矩阵的稀疏矩阵函数调用格式
$a = [a_1, a_2, \cdots, a_m, n];$
$b = [b_1, b_2, \cdots, b_m, n];$
$w = [w_1, w_2, \cdots, 2_m, 0];$
R = sparse(a, b, w);
$其中，a_i 为起始节点向量，b_i为终止节点向量，w_i为边权值向量，这里i = 1, 2, \cdots, m$
各个向量最后的一个值使得关联矩阵为方阵

% 建立关联矩阵并显示图
ab = [1 1 2 2 3 3 4 4 4 4 5 6 6 7 8];
bb = [2 3 5 4 4 6 5 7 8 6 7 8 9 9 9];
w = [1 2 12 6 3 4 4 15 7 2 7 7 15 3 10];
R = sparse(ab, bb, w);
R(9, 9) = 0;
h = view(biograph(R, [], 'ShowWeights', 'on'))
% 搜索最短路径
[d, p] = graphshortestpath(R, 1, 9)
set(h.Nodes(p), 'Color', [1, 0.4 0.4])
edges = getedgesbynodeid(h, get(h.Nodes(p), 'ID'));
set(edges, 'LineColor', [1 0 0])

无向图的路径最优搜索
无向图的关联矩阵构造方法
先按照有向图的方式构造关联矩阵R
无向图的关联矩阵为 $R_1 = R + R^T$
如果无向图中，某些边是有向的，则手工修改该矩阵
若由节点i到节点j与由节点j到i的边权值是不同的，用手工方法重新定义和修改