数学规划是运筹学的一个分支，研究在给定的条件下（约束条件），如何按照某一标准（目标函数）寻求计划的最佳方案。也就是求目标函数在一定约束条件下的极值的问题。不同类型的规划问题需要采用不同的建模方法和求解算法，因此在解决具体问题时需要根据问题特点选择合适的方法进行建模和求解。

数学规划问题通常可以分为以下几类：

  • 线性规划（Linear Programming, LP）：目标函数和约束条件均为线性的规划问题。例如资源分配、生产计划、运输优化等。

  • 非线性规划（Nonlinear Programming, NLP）：目标函数和/或约束条件包含非线性部分的规划问题。求解相对复杂，涉及到局部最优解和全局最优解的搜索。

  • 整数规划（Integer Programming，IP）：在线性规划的基础上，自变量需要取整数值的优化问题。

  • 二次规划（Quadratic Programming，QP）：目标函数为二次函数，约束条件为线性函数的优化问题。

我们主要讨论泛用性较高的线性规划和非线性规划。

线性规划

要解决线性规划问题，首先需要定义决策变量、目标函数和约束条件。然后，可以使用MATLAB中的优化函数对问题进行求解。

举一个例子,要求解的线性规划问题的目标是最大化目标函数 f = -5a - 4b - 6c，其中a、b、c是决策变量。

约束条件为：

a - b + c ≤ 20
3a + 2b + 4c ≤ 42
3a + 2b ≤ 30
a、b、c ≥ 0

为了解决这个问题，我们将调用 linprog 函数进行线性规划求解，代码是：

1
[x, fval] = linprog(c, A, b, Aeq, beq, lb, ub)

其中的参数含义如下：

  • c：目标函数的系数向量。它是一个列向量，表示目标函数中每个决策变量的系数。

  • A：不等式约束条件左侧的系数矩阵。它是一个矩阵，每一行对应一个不等式约束条件，并且每一列对应一个决策变量的系数。

  • b：不等式约束条件右侧的常数向量。它是一个列向量，表示不等式约束条件的右侧值。

  • Aeq：等式约束条件左侧的系数矩阵。它是一个矩阵，每一行对应一个等式约束条件，并且每一列对应一个决策变量的系数。

  • beq：等式约束条件右侧的常数向量。它是一个列向量，表示等式约束条件的右侧值。

  • lb：决策变量的下界向量。它是一个列向量，表示每个决策变量的下限值。如果某个决策变量没有下限，则对应的值设为 -Inf，如果不写就是没有。

  • ub：决策变量的上界向量。它是一个列向量，表示每个决策变量的上限值。如果某个决策变量没有上限，则对应的值设为 Inf,如果不写就是没有。

要使用这个代码，我们就需要把问题表示成线性规划问题的标准型。在MATLAB中，线性规划问题的标准型指的是形如下面形式的线性规划问题：

1
2
3
4
5
minimize c^T * x
subject to:
    A * x <= b
    Aeq * x = beq
    lb <= x <= ub

这里都是“小于”形式，包括求解也是求最小值。如果有“大于”的形式出现，则等式两边同时加一个负号就可以了。

如果约束条件是大于，在约束条件里加上负号。

如果整个问题求的是最大值，我们先在c（目标函数的系数向量）里把目标函数的系数全部变号，再把最后求解的值变号即可。也就是说，**max f(x) = min -f(x)**。

解决上述例子的代码是：

1
2
3
4
5
6
7
8
9
10
c = [5 4 6]';  % 加单引号表示转置，这里需要变号

A = [1 -1 1;
     3 2 4;
     3 2 0];
b = [20 42 30]';   
lb = [0 0 0]'; 
[x fval] = linprog(c, A, b, [], [], lb) 
fval = -fval % 要取负号，原来是求最大值，我们添加负号变成了max f(x) = min -f(x)

例子随手改的，答案全是0。

部分特殊情况：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
%% 多个解的情况
% 例如 ： min z = x1 + x2   s.t.  x1 + x2 >= 10
c = [1 1]';   
A = [-1 -1];
b = -10;
[x fval] = linprog(c, A, b)   % Aeq, beq, lb和ub我们都没写，意味着没有等式约束和上下界约束
% x有多个解时，Matlab会给我们返回其中的一个解

%% 不存在解的情况
% 例如 ： min z = x1 + x2   s.t.  x1 + x2 = 10 、 x1 + 2*x2 <= 8、 x1 >=0 ，x2 >=0 
c = [1 1]'; 
A = [1 2];
b = 8;
Aeq = [1 1];
beq = 10;
lb = [0 0]';
[x fval] = linprog(c, A, b, Aeq, beq, lb)  % Linprog stopped because no point satisfies the constraints.（没有任何一个点满足约束条件）

非线性规划

很显而易见的是，非线性规划是指目标函数或约束条件中包含非线性项的优化问题。也就是说，在非线性规划中，目标函数和约束条件可以是非线性的，即涉及到变量的非线性函数、乘积、幂等运算或其他非线性形式。

举一个例子,要求解的非线性规划问题的目标是最大化目标函数f(x) = x1^2 + x2^2 - x1*x2 - 2x1 - 5x2。

约束条件为：

-(x1-1)^2 + x2 >= 0
2x1 - 3x2 + 6 >= 0

这两个不等式分别表示了在定义域内 x1 和 x2 的取值范围限制。

同样的，为了求解非线性规划问题，我们也需要把问题表示成非线性规划问题的标准型:

1
2
3
4
5
minimize c^T * x
subject to:
    A * x <= b , Aeq * x = beq; (线性不等式约束)
    c(x) <= 0 , ceq(x) = 0; （非线性不等式约束）
    lb <= x <= ub;（上界，下界）

对于非线性规划，我们使用fmincon 函数，标准格式为：

1
[x, fval, exitflag, output] = fmincon(fun, x0, A, b, Aeq, beq, lb, ub, nonlcon, options)

其中的参数含义如下：

  • fun：目标函数，该函数接受一个自变量向量作为输入，并返回一个标量作为输出。

  • x0：自变量的初始值，是一个向量。与线性规划相比，再非线性规划问题里，对它的选取非常主要，因为非线性规划的算法求解是一个局部最优值。

  • A 和 b：不等式约束的系数矩阵和右侧系数。A 是一个 m×n 的矩阵，b 是一个 m×1 的向量，表示 m 个不等式约束的系数和右侧限制。

  • Aeq 和 beq：等式约束的系数矩阵和右侧系数。Aeq 是一个 p×n 的矩阵，beq 是一个 p×1 的向量，表示 p 个等式约束的系数和右侧限制。

  • lb 和 ub：自变量的下界和上界，分别是 n×1 的向量。若没有下界或上界，则使用 -inf 或 inf 表示。

  • nonlcon：非线性约束函数，该函数接受一个自变量向量作为输入，并返回一个向量作为输出，表示非线性约束的结果。

  • options：优化选项，可以使用 optimoptions 函数创建一个选项对象，也可以使用 optimset 函数创建一个选项结构体。

需要注意，这里的fun和nonlcon是需要自己写成函数形式的。

对于fun,可以任意取名，到时候会被fmincon函数调用, 它实际上就是目标函数,输入值x实际上就是决策变量，由x1和x2组成的向量，代码如下：

1
2
3
function f = fun1(x)
    f = -x(1)^2-x(2)^2 +x(1)*x(2)+2*x(1)+5*x(2) ; 
end

这里的f实际上就是目标函数，函数的返回值也是f。fun1是函数名称，保存的m文件和函数名称得一致，也要为fun1.m。

对于nonlcon，和fun类似，输入值x实际上就是决策变量，由x1和x2组成的一个向量。返回值有两个，一个是非线性不等式约束c，一个是非线性等式约束ceq。nonlfun1是函数名称，到时候会被fmincon函数调用。保存的m文件和函数名称得一致，也要为nonlfun1.m

1
2
3
4
function [c,ceq] = nonlfun1(x)
   c = [(x(1)-1)^2-x(2)];   %x1 写成 x(1)
   ceq = [];  % 不存在非线性等式约束，所以用[]表示
end

现在还有一个问题，就是初始值的选择。初始值指的是优化算法开始搜索的起点，这个起点可以是一个参数向量、一组变量的取值，或者是一个解空间中的初始点。优化算法会根据这个初始值开始搜索最优解或局部最优解。如果选择的初始值不好，可能会导致优化算法无法在可接受的时间内找到全局最优解。在更极端的情况下，选择不当的初始值可能导致优化算法陷入无限循环或者发散。

在这里，我们使用蒙特卡洛算法处理。完整代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
clc,clear;
n=10000000; %生成的随机数组数
x1=unifrnd(-100,100,n,1);  % 生成在[-100,100]之间均匀分布的随机数组成的n行1列的向量构成x1
x2=unifrnd(-100,100,n,1);  % 生成在[-100,100]之间均匀分布的随机数组成的n行1列的向量构成x2，这两步具体情况具体分析
fmin=+inf; % 初始化函数f的最小值为正无穷,后续再更新
for i=1:n
    x = [x1(i), x2(i)];  %构造x向量, 这里千万别写成了：x =[x1, x2]
    if ((x(1)-1)^2-x(2)<=0)  & (-2*x(1)+3*x(2)-6 <= 0)     % 判断是否满足约束条件
        result = -x(1)^2-x(2)^2 +x(1)*x(2)+2*x(1)+5*x(2) ;  % 如果满足条件，就计算函数值
        if  result  < fmin  
            fmin = result;  
            x0 = x;  % 更新
        end
    end
end
disp('蒙特卡罗选取的初始值为：'); disp(x0)

A = [-2 3]; b = 6;
[x,fval] = fmincon(@fun1,x0,A,b,[],[],[],[],@nonlfun1)
fval = -fval