数学建模求解彩票选择问题

下面来介绍一下如何在MATLAB进行数学建模，来实现对彩票问题的组合选择求解，具体如下：

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

1、在MATLAB主界面的编辑器中分别写入下列代码：

主程序：

% 本程序包含的m文件为：
%     main.m：主程序
%     cpiao.m：目标函数
%     calculate_probability.m：计算各奖项获奖概率
%     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

% 对于K2、K3、K4型的情况
% n选m或(m+1)，n的选择范围在29到60，m的选择范围为5到7
% 故有 (60-29+1)*(7-5+1)=96种 取法，
% 依题意，K2、K3、K4都有这96种取法，也即第三维上为3
% 故有下面的变量声明：
p_all=zeros(7,96,3);
rx_all = zeros(7,96,3);
fval_all= zeros(1,96,3);
flag_all = zeros(1,96,3);
for m=5:7
    for n=29:60
        for i=1:3
            % 根据i的值判断属于（K2、K3、K4）中哪一型
            % (i=1是K2;i=2是K3;i=3是K4)，
            % 并根据n、m生成各奖项概率
            % p_temp=eval(sprintf('comb_k%d(m,n)',i+1));
			p_temp = calculate_probability(m,n,i+1);
            p_all(:,(m-5).*32+(n-28),i) = p_temp;
            % K2、K3可合并处理（奖项数目一样）
            if (i ~= 3)
                Aeq=[1,1,1,0,0,0,0];
                beq=1;
                a_lb=[10,4,3,4,2,2];
                b_ub=[233,54,17,20,10,10];
                A=[ 0,0,0,-1,a_lb(4),0,0;
                    0,0,0,1,-b_ub(4),0,0;
                    0,0,0,0,-1,a_lb(5),0;
                    0,0,0,0,1,-b_ub(5),0];
                %由于x(7)可能为零，故不在这里对x(6)/x(7)进行上下限限制，
                % 而在非线性约束nonlcon中进行
                %    0,0,0,0,0,-1,a_lb(6);
                %    0,0,0,0,0,1,-b_ub(6)