数学建模（一）规划问题

1.线性规划

在一组线性的约束条件的限制下，求一线性目标函数最大或最小的问题。
直接上函数: linprog(c,A,b)

形如[x,fval]=linprog(c,A,b,Aeq,beq,LB,UB,X0,OPTIONS)

c=[2;3;-5];
a=[-2,5,-1;1,3,1];
b=[-10;12];
aeq=[1,1,1];
beq=7;
x=linprog(-c,a,b,aeq,beq,zeros(3,1))

Optimal solution found.


x =

    6.4286
    0.5714
         0

>> value=c'*x

value =

   14.5714

对应着来写就可以了，只有一个点要说就是函数标准默认求的是最小值，如果求最大值只需把目标函数取负，like代码里的-c。
Lingo更简单一点：

max=2*x1+3*x2-5*x3;
x1+x2+x3=7;
2*x1-5*x2+x3>=10;
x1+3*x2+x3<=12;
x1>=0;
x2>=0;
x3>=0;

例：不例了，还是套公式，重点是怎么从实际问题抽象出模型，且行且珍惜吧。

2.整数规划

定义：规划中的变量（部分或全部）限制为整数时。
LINGO之。一个例子：

model:
sets:
row/1..4/:b;
col/1..5/:c1,c2,x;
link(row,col):a;
endsets
data:
c1=1,1,3,4,2;
c2=-8,-2,-3,-1,-2;
a=1 1 1 1 1
1 2 2 1 6 
2 1 6 0 0
0 0 1 1 5;
b=400,800,200,200;
enddata
max=@sum(col:c1*x^2+c2*x);
@for(row(i):@sum(col(j):a(i,j)*x(j))<b(i));
@for(col:@gin(x));
@for(col:@bnd(0,x,99));
end 

@gin是整数限定。但我不明白为啥写成这么复杂，简化之：

max=x1^2+x2^2+3*x3^2+4*x4^2+2*x5^2-8*x1-2*x2-3*x3-x4-2*x5;
x1<=99;
x2<=99;
x3<=99;
x4<=99;
x5<=99;
x1>=0;
x2>=0;
x3>=0;
x4>=0;
x5>=0;
x1+x2+x3+x4+x5<=400;
x1+2*x2+2*x3+x4+6*x5<=800;
2*x1+x2+6*x3<=200;
x3+x4+5*x5<=200;
@gin(x1);
@gin(x2);
@gin(x3);
@gin(x4);
@gin(x5);

好像并没有简单。

3.非线性规划

目标函数或约束条件中包含非线性函数。
函数：X=FMINCON(FUN,X0,A,B,Aeq,Beq,LB,UB,NONLCON,OPTIONS)。

例：

首先需要编写一个目标函数文件

function f=fun1(x);
f=sum(x.^2)+8;

然后需要编写一个约束函数文件

function [g,h]=fun2(x);
g=[-x(1)^2+x(2)-x(3)^2
x(1)+x(2)^2+x(3)^3-20]; %非线性不等式约束
h=[-x(1)-x(2)^2+2
x(2)+2*x(3)^2-3];

然后即可执行

options=optimset('largescale','off');
>> [x,y]=fmincon('fun1',rand(3,1),[],[],[],[],zeros(3,1),[], ...
'fun2', options)

Local minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the default value of the optimality tolerance,
and constraints are satisfied to within the default value of the constraint tolerance.

<stopping criteria details>


x =

    0.5522
    1.2033
    0.9478


y =

   10.6511

可以自己定义一个函数，感觉方便多了。
lingo之：

min=x1^2+x2^2+x3^2+8;
x1^2-x2+x3^2>=0;
x1+x2^2+x3^3<=20;
-x1-x2^2+2=0;
x2+2*x3^2=3;
x1>=0;
x2>=0;
x3>=0;

答案一样的。
例:CUMCM 2002B 彩票中的数学
因为这个题就是对多种中奖方案的目标函数求最大值然后选择方案的问题，是一个比较典型的(非线性)规划型问题。首先给出建立的规划模型：

代码只粘主函数的k1部分了，注释也挺清楚的，需要注意的就是为了避免求得的是局部最优解，可以多次设初始值。

% 功能说明：
%     根据参考答案中的模型，本程序分别对K1、K2、K3、K4型彩票进行求解，并对
%     n、m的各种组合进行循环。求解时，首先计算当前n、m的各奖项获奖概率，然后
%     随机生成多个初始值，调用fmincon函数寻找目标函数的最小值（原目标函数要求
%     极大，但fmincon是寻找极小，故令原目标函数乘以(-1)，寻找新目标函数的极小值），
%     最后比较各种类型彩票的求解结果，输出对应最大的原目标函数的解。
% 本程序包含的m文件为：
%     main.m：主程序
%     cpiao.m：目标函数
%     calculate_probability.m：计算各奖项获奖概率也就是目标函数里的p
%     nonlcon.m：非线性约束
% 使用说明：
%    执行main.m

clc 
clear 

% 为避免陷入局部最优，需要以随机的初值进行多次尝试，
% 该变量为对每个m/n组合生成随机初值的数目，越大则找到
% 全局最优的概率越大，但程序运行的时间也越长，
% 请根据电脑情况自行设置
nums_test_of_initial_value = 20; 

global v
v = 630589;		% 求解v为630589的收入水平情况
DEBUG = 0;
rand('state',sum(100*clock))	% 初始化随机数生成器
format long g

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 求解开始
% 对于K1型
p_k1 = [2e-7;8e-7;1.8e-5;2.61e-4;3.42e-3;4.1995e-2];
% 6个奖项6个变量
Aeq=[1,1,1,0,0,0];
beq=1;
a_lb=[10,4,3,4,2];
b_ub=[233,54,17,20,10];
A= [0,0,0,-1,a_lb(4),0;
    0,0,0,1,-b_ub(4),0;
    0,0,0,0,-1,a_lb(5);
    0,0,0,0,1,-b_ub(5)];
b= [0;0;0;0];
lb=[0.5;0;0;0;0;0];
ub=[0.8;1;1;inf;inf;inf];
p_test = p_k1;
rx0_tmp = zeros(6,1);%初始值
rx_meta_result = zeros(6,1);
fval_meta_result = inf;
flag_meta_result = nan; %用以判断有没有得到过可行解
if DEBUG == 1
    output_meta_result = [];
end
for j = 1:nums_test_of_initial_value
    %随机生成多个初始值rx0_tmp，以避免局部最优
    rx0_tmp(1) = rand*(0.8-0.5) + 0.5;
    rx0_tmp(2) = rand*(1-rx0_tmp(1));
    rx0_tmp(3) = 1 - rx0_tmp(1) - rx0_tmp(2);
    rx0_tmp(4) = rand*1000;
    rx0_tmp(5) = rand*100;
    rx0_tmp(6) = rand*50;
    % 寻优
    [rx_tmp,fval_tmp,flag_tmp,output_tmp]= ...
            fmincon('cpiao',rx0_tmp,A,b,...
                    Aeq,beq,lb,ub,'nonlcon',[],1,p_test,a_lb,b_ub);
	% 上式倒数第四个参数是为了区分彩票的类型(K1/K2/K3/K4)
	% 最后三个是函数cpiao和nonlcon计算中可能要用到的量。
    if (flag_tmp == 1) && (fval_meta_result > fval_tmp)
        fval_meta_result = fval_tmp;
        rx_meta_result = rx_tmp;
        flag_meta_result = 1;
        if DEBUG == 1
            output_meta_result = output_tmp;
        end
    end
end
% 把求得的最好结果保存下来
if ~isnan(flag_meta_result)
    rx_k1 = rx_meta_result;
    fval_k1 = fval_meta_result;
    flag_k1 = flag_meta_result;
    if DEBUG == 1
        output = output_meta_result;
    end
else
    if DEBUG == 1
        rx_k1 = rx_tmp;
        fval_k1 = fval_tmp;
        flag_k1 = flag_tmp;
        output = output_tmp;
    end
end

这个在设变量的时候不要那么死板的照着公式写，约束函数看着挺多有10个，其实写到约束函数文件里的只有式（6），其他的式子等于都隐含在目标函数或是初始值等地方了。

4.动态规划

把多阶段过程转化为一系列单阶段问题，逐个求解。主要用于求解以时间划分阶段的动态过程的优化问题。但是，一些与时间无关的静态规划（如线性规划和非线性规划），只要人为引进时间因素，把它视为多阶段决策过程，亦可以用动态规划求解。
动态规划可以求解的问题类型很多，这里依然举一个简单的例子：最短路径问题。下图是一个线路网， 连线上的数字表示两点之间的距离（或费用）。试寻求一条由 A到G 距离最短（或费用最省）的路线。

求解步骤:
（ i）将过程划分成恰当的阶段。
（ ii）正确选择状态变量 xk ，使它既能描述过程的状态，又满足无后效性，同时确定允许状态集合 Xk 。
（ iii）选择决策变量uk ，确定允许决策集合Uk (xk ) 。
（ iv）写出状态转移方程。
（ v）确定阶段指标 vk (xk ,uk ) 及指标函数Vkn 的形式（阶段指标之和，阶段指标之积，阶段指标之极大或极小等）。
（ vi）写出基本方程即最优值函数满足的递归方程，以及端点条件。
Lingo解之：

model:
Title Dynamic Programming;
sets:
vertex/A,B1,B2,C1,C2,C3,C4,D1,D2,D3,E1,E2,E3,F1,F2,G/:L;
road(vertex,vertex)/A B1,A B2,B1 C1,B1 C2,B1 C3,B2 C2,B2 C3,B2 C4,
C1 D1,C1 D2,C2 D1,C2 D2,C3 D2,C3 D3,C4 D2,C4 D3,
D1 E1,D1 E2,D2 E2,D2 E3,D3 E2,D3 E3,
E1 F1,E1 F2,E2 F1,E2 F2,E3 F1,E3 F2,F1 G,F2 G/:D;
endsets
data:
D=5 3 1 3 6 8 7 6
6 8 3 5 3 3 8 4
2 2 1 2 3 3
3 5 5 2 6 6 4 3;
L=0,,,,,,,,,,,,,,,;       ！15个点，大概代表15个城市
enddata
@for(vertex(i)|i#GT#1:L(i)=@min(road(j,i):L(j)+D(j,i)));

我看到后面启发式算法时举的例子是TSP旅行商问题，跟这个例子很像，如果用lingo就能这么简单的解出来，那为何还要费劲用什么模拟退火呢？？？
这个问题留到后面比较一下。
综上规划问题就统一用lingo吧，smile。