Maltab在数学建模中的应用（第二版）——读书笔记下

最新推荐文章于 2020-10-26 14:49:46 发布
最新推荐文章于 2020-10-26 14:49:46 发布
阅读量1k 收藏 20
点赞数 4
分类专栏： matlab 文章标签： 数学建模 matlab
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/qq_44809326/article/details/108338883
版权

Maltab在数学建模中的应用（第二版）——读书笔记下

数学建模问题可分为三大类：预测类问题、连续性优化问题、离散型优化问题。

预测类、连续性优化问题要求有很强的数学基础，但是求解方法固定。

离散型优化问题难在求解，所以注重各种优化算法用matlab自主编程求解规划问题。

  1. 预测类问题
    使用数学方法：拟合、回归、插值、神经网络、灰色预测、小波分析
    CUMCM：2003A、2005A、2007A

  2. 连续型优化问题
    使用数学方法：拟合、回归、插值、微分、求极值
    CUMCM：2002A、2008A、2009A

  3. 离散型优化问题

    使用数学方法：目标规划

    CUMCM：2002B、2003B、2004A、2004B、2005B、2006A、2007B、2009B

1彩票中的数学2002B

文档：2002B.note
链接：http://note.youdao.com/noteshare?id=39677d3761846916df202611a3c6e104&sub=92EE275E5AC64029B7D7683EEA177130

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


%1.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)];
                %b=[0;0;0;0;0;0];
                b=[0;0;0;0];
                lb=[0.5;0;0;0;0;0;0];
                ub=[0.8;1;1;inf;inf;inf;inf];
                p_test = p_temp;
                %随机生成多个初始值rx0_tmp，以避免局部最优
                rx0_tmp = zeros(7,1);
                rx_meta_result = zeros(7,1);
                fval_meta_result = inf;
                flag_meta_result = nan; %用以判断有没有得到过可行解
                for j = 1:nums_test_of_initial_value
                    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;
                    rx0_tmp(7) = rand*10;
                    [rx_tmp,fval_tmp,flag_tmp]= ...
                            fmincon('cpiao',rx0_tmp,...
                                    A,b,Aeq,beq,lb,ub,...
										'nonlcon',[],i+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;
                    end
                end
                % 把求得的最好结果保存下来
                rx_all(:,(m-5).*32+(n-28),i) = rx_meta_result;
                fval_all(1,(m-5).*32+(n-28),i) = fval_meta_result;
                flag_all(1,(m-5).*32+(n-28),i) = flag_meta_result;
            else
            % i==3，相应于K4型
                % 对于K4型，因只设到五等奖，故仅5个变量了
                Aeq=[1,1,1,0,0];
                beq=1;
                a_lb=[10,4,3,4];
                b_ub=[233,54,17,20];
                A=[ 0,0,0,-1,a_lb(4);
                    0,0,0,1,-b_ub(4)];
                b=[0;0];
                lb=[0.5;0;0;0;0];
                ub=[0.8;1;1;inf;inf];
                p_test = p_temp;
                %随机生成多个初始值rx0_tmp，以避免局部最优
                rx0_tmp = zeros(5,1);
                rx_meta_result = zeros(5,1);
                fval_meta_result = inf;
                flag_meta_result = nan; %用以判断有没有得到过可行解
                for j = 1:nums_test_of_initial_value
                    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;
                    [rx_tmp,fval_tmp,flag_tmp]= ...
                            fmincon('cpiao',rx0_tmp,A,b,...
                                    Aeq,beq,lb,ub,'nonlcon',...
										[],4,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;
                    end
                end
                % 把求得的最好结果保存下来
                rx_all(:,(m-5).*32+(n-28),i) = [rx_meta_result;0;0];
                fval_all(1,(m-5).*32+(n-28),i) = fval_meta_result;
                flag_all(1,(m-5).*32+(n-28),i) = flag_meta_result;
                
            end
        end
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 寻优结束，进行结果处理
% 在所有（K1、K2、K3、K4）求解结果中找目标函数最小的

% 判断（K1、K2、K3、K4）中哪一种的目标函数最小
ind_tmp = (flag_all >= 0);
if sum(sum(sum(ind_tmp))) ~= 0
    % K2、K3、K4的求解中找到了可行解（或最优解）
    val_tmp = fval_all.*ind_tmp;
    [val_tmp2,ind_tmp2] = min(val_tmp);
    [val_min,ind_tmp3] = min(val_tmp2);
    if (flag_k1 < 0)
        % K1的求解中没找到可行解
        signal = 1;     % 标志变量
    else
        if val_min < fval_k1
            signal = 1;
        elseif val_min > fval_k1
            signal = 2;
        else
            signal = 3;
        end
    end
    ;
else
    % K2、K3、K4的求解中没有找到了可行解
    if (flag_k1 < 0)
        % K1的求解中没找到可行解
        disp('(K1、K2、K3、K4)所有的求解中')
        disp('一个可行解都没有找到')
        disp('(还并不意味着完全没有可行解，')
        disp('也许是初值点选的不好因此没有找到）')
        break;
    else
        signal = 2;
    end
end

if (signal == 1)
    ind_tmp4 = ind_tmp2(ind_tmp3);
    rx_result = rx_all(:,ind_tmp4,ind_tmp3);
    fval_result = fval_all(:,ind_tmp4,ind_tmp3);
    fval_result =-fval_result;
    n = (ind_tmp4 - floor(ind_tmp4 / 32) * 32) + 28;
    m = floor(ind_tmp4 / 32) + 5;
    p_tmp = p_all(:,ind_tmp4,ind_tmp3);
elseif signal == 2
    rx_result = rx_k1;
    fval_result =-fval_k1;
    p_tmp = p_k1;
else    %signal == 3
    ind_tmp4 = ind_tmp2(ind_tmp3);
    rx_result = rx_all(:,ind_tmp4,ind_tmp3);
    fval_result = fval_all(:,ind_tmp4,ind_tmp3);
    fval_result =-fval_result;
    n = (ind_tmp4 - floor(ind_tmp4 / 32) * 32) + 28;
    m = floor(ind_tmp4 / 32) + 5;
end

% 输出计算结果
if signal == 1      % 最优解在K2、K3、K4中时
    if ind_tmp3 == 1
        disp(sprintf('最优解为：K2型，%d选%d',n,m));
    elseif ind_tmp3 == 2
        disp(sprintf('最优解为：K3型，%d选%d+1',n,m));
    elseif ind_tmp3 == 3
        disp(sprintf('最优解为：K4型，%d选%d无特别号',n,m));
    end
elseif signal == 2  % 最优解在K1中时
    disp(sprintf('最优解为：K1型，10选6+1'));
else                % K1的解和K2、K3、K4的解重合时
    if ind_tmp3 == 1
        disp(sprintf('10选6+1和K2型%d选%d同为最优解',n,m));
    elseif ind_tmp3 == 2
        disp(sprintf('10选6+1和K3型%d选%d+1同为最优解',n,m));
    elseif ind_tmp3 == 3
        disp(sprintf('10选6+1和K4型%d选%d无特别号同为最优解',n,m));
    end
end
disp('对应的目标函数值为：')
disp(fval_result)
if signal ~= 3
    disp('最终求解变量值为：')
    disp(rx_result)
    disp('各奖项的金额是：')
    x=zeros(3,1);
    x(1) = (1-p_tmp(4).*rx_result(4)-...
              p_tmp(5).*rx_result(5)-...
              p_tmp(6).*rx_result(6)-...
              p_tmp(7).*rx_result(7)).*rx_result(1)./p_tmp(1);
    x(2) = (1-p_tmp(4).*rx_result(4)-...
              p_tmp(5).*rx_result(5)-...
              p_tmp(6).*rx_result(6)-...
              p_tmp(7).*rx_result(7)).*rx_result(2)./p_tmp(2);
    x(3) = (1-p_tmp(4).*rx_result(4)-...
              p_tmp(5).*rx_result(5)-...
              p_tmp(6).*rx_result(6)-...
              p_tmp(7).*rx_result(7)).*rx_result(3)./p_tmp(3);
    rx_money=[x;rx_result(4:7)];
    disp(rx_money)
else
    disp('最终求解变量值为：')
    disp('10选6+1时')
    disp(rx_k1)
    disp('K%d型时',ind_tmp3+1)
    disp(rx_result)
end




%2.calculate_probability.m
function p_temp_sub = calculate_probability(m,n,type);
%n选m时各奖项的获奖概率
% type代表当前求解的彩票类型
% K1型是固定概率，无需在这里计算

if (type == 2)
	% K2型的
	p_temp_sub=zeros(7,1);
	p_temp_sub(1)=1/mmmcomb(n,m);
	p_temp_sub(2)=mmmcomb(m,m-1)./mmmcomb(n,m);
	p_temp_sub(3)=mmmcomb(m,m-1).*mmmcomb(n-m-1,1)./mmmcomb(n,m);
	p_temp_sub(4)=mmmcomb(m,m-2).*mmmcomb(n-m-1,1)./mmmcomb(n,m);
	p_temp_sub(5)=mmmcomb(m,m-2).*mmmcomb(n-m-1,2)./mmmcomb(n,m);
	p_temp_sub(6)=mmmcomb(m,m-3).*mmmcomb(n-m-1,2)./mmmcomb(n,m);
	p_temp_sub(7)=mmmcomb(m,m-3).*mmmcomb(n-m-1,3)./mmmcomb(n,m);
elseif (type == 3)
	% K3型的
	p_temp_sub=zeros(7,1);
	p_temp_sub(1)=1./mmmcomb(n,m+1);
	p_temp_sub(2)=mmmcomb(n-m-1,1)./mmmcomb(n,m+1);
	p_temp_sub(3)=mmmcomb(m,m-1).*mmmcomb(n-m-1,1)./mmmcomb(n,m+1);
	p_temp_sub(4)=mmmcomb(m,m-1).*mmmcomb(n-m-1,2)./mmmcomb(n,m+1);
	p_temp_sub(5)=mmmcomb(m,m-2).*mmmcomb(n-m-1,2)./mmmcomb(n,m+1);
	p_temp_sub(6)=mmmcomb(m,m-2).*mmmcomb(n-m-1,3)./mmmcomb(n,m+1);
	p_temp_sub(7)=mmmcomb(m,m-3).*mmmcomb(n-m-1,3)./mmmcomb(n,m+1);
elseif (type == 4)
	% K4型的
	p_temp_sub=zeros(7,1);
	p_temp_sub(1)=1./mmmcomb(n,m);
	p_temp_sub(2)=mmmcomb(m,m-1).*mmmcomb(n-m,1)./mmmcomb(n,m);
	p_temp_sub(3)=mmmcomb(m,m-2).*mmmcomb(n-m,2)./mmmcomb(n,m);
	p_temp_sub(4)=mmmcomb(m,m-3).*mmmcomb(n-m,3)./mmmcomb(n,m);
	p_temp_sub(5)=mmmcomb(m,m-4).*mmmcomb(n-m,4)./mmmcomb(n,m);
	p_temp_sub(6)=0;
	p_temp_sub(7)=0;
else
	error('Error in calculate_probability!');
end


function combi = mmmcomb(n,m)
%求从n个数中取出m个数的组合数
if (isscalar(n)) && (isscalar(m)) &&...
		(isreal(n)) && (isreal(m)) && (n>=m) && (m>0)
    combi=factorial(n)./factorial(m)./factorial(n-m);
else
    error('A mistake occurs when calculating combinations.')
end




%3.cpiao.m
function f=cpiao(rx,type,p_test,a_lb,b_ub)
% type代表当前求解的彩票类型
% p_test是当前各奖项概率（当前彩票类型，当前m、n）
% a_lb,b_ub在本函数中用不到
% 注意：这里的目标函数是原目标函数乘以 -1；

global v
if (type == 1)
	% 这是K1型的
	rx_last3=rx(4:6);
	p_last3=p_test(4:6);
	p_last3=p_last3';
	sum_last3=p_last3*rx_last3;
	f= ...
	-(  p_test(1).*(1-exp(-(((1-sum_last3).*rx(1))./p_test(1)./v).^2))+...
		p_test(2).*(1-exp(-(((1-sum_last3).*rx(2))./p_test(2)./v).^2))+...
		p_test(3).*(1-exp(-(((1-sum_last3).*rx(3))./p_test(3)./v).^2))+...
		p_test(4).*(1-exp(-(rx(4)./v).^2))+...
		p_test(5).*(1-exp(-(rx(5)./v).^2))+...
		p_test(6).*(1-exp(-(rx(6)./v).^2)) ...
	 );
elseif (type == 2) || (type == 3)
	% 这是K2和K3的（K2、K3可合并处理，因奖项数目一样）
	rx_last4=rx(4:7);
	p_last4=p_test(4:7);
	p_last4=p_last4';
	sum_last4=p_last4*rx_last4;
	f= ...
	-(  p_test(1).*(1-exp(-(((1-sum_last4).*rx(1))./p_test(1)./v).^2))+...
		p_test(2).*(1-exp(-(((1-sum_last4).*rx(2))./p_test(2)./v).^2))+...
		p_test(3).*(1-exp(-(((1-sum_last4).*rx(3))./p_test(3)./v).^2))+...
		p_test(4).*(1-exp(-(rx(4)./v).^2))+...
		p_test(5).*(1-exp(-(rx(5)./v).^2))+...
		p_test(6).*(1-exp(-(rx(6)./v).^2))+...
		p_test(7).*(1-exp(-(rx(7)./v).^2)) ...
	 );
elseif (type == 4)
	% K4型
	rx_last2=rx(4:5);
	p_last2=p_test(4:5);
	p_last2=p_last2';
	sum_last2=p_last2 * rx_last2;
	f= ...
	-(  p_test(1).*(1-exp(-(((1-sum_last2).*rx(1))./p_test(1)./v).^2))+...
		p_test(2).*(1-exp(-(((1-sum_last2).*rx(2))./p_test(2)./v).^2))+...
		p_test(3).*(1-exp(-(((1-sum_last2).*rx(3))./p_test(3)./v).^2))+...
		p_test(4).*(1-exp(-(rx(4)./v).^2))+...
		p_test(5).*(1-exp(-(rx(5)./v).^2)) ...
	 );
else
	error('Error in function cpiao!')
end

%4.nonlcon.m
function [c,ceq]=nonlcon(rx,type,p_test,a_lb,b_ub)
% type代表当前求解的彩票类型
% p_test是当前各奖项概率（当前彩票类型，当前m、n）
% a_lb,b_ub是相邻两个奖项奖金之比的下限和上限

if (type == 1)
	% 这是K1型的
	c(1)= 6e5- ...
			 (1-p_test(4).*rx(4)-p_test(5).*rx(5)-...
				p_test(6).*rx(6)).*rx(1)./p_test(1);
	c(2)=-5e6+ ...
			 (1-p_test(4).*rx(4)-p_test(5).*rx(5)-...
				p_test(6).*rx(6)).*rx(1)./p_test(1);
	c(3)= a_lb(1).* ...
			 (1-p_test(4).*rx(4)-p_test(5).*rx(5)-...
				p_test(6).*rx(6)).*rx(2)./p_test(2)-...
			 (1-p_test(4).*rx(4)-p_test(5).*rx(5)-...
				p_test(6).*rx(6)).*rx(1)./p_test(1);
	c(4)= a_lb(2).* ...
			 (1-p_test(4).*rx(4)-p_test(5).*rx(5)-...
				p_test(6).*rx(6)).*rx(3)./p_test(3)-...
			 (1-p_test(4).*rx(4)-p_test(5).*rx(5)-...
				p_test(6).*rx(6)).*rx(2)./p_test(2);
		
	c(5)= a_lb(3).* rx(4)-...
			 (1-p_test(4).*rx(4)-p_test(5).*rx(5)-...
				p_test(6).*rx(6)).*rx(3)./p_test(3);

	c(6)=-b_ub(1).* ...
			 (1-p_test(4).*rx(4)-p_test(5).*rx(5)-...
				p_test(6).*rx(6)).*rx(2)./p_test(2)+...
			 (1-p_test(4).*rx(4)-p_test(5).*rx(5)-...
				p_test(6).*rx(6)).*rx(1)./p_test(1);
		
	c(7)=-b_ub(2).* ...
			 (1-p_test(4).*rx(4)-p_test(5).*rx(5)-...
				p_test(6).*rx(6)).*rx(3)./p_test(3)+...
			 (1-p_test(4).*rx(4)-p_test(5).*rx(5)-...
				p_test(6).*rx(6)).*rx(2)./p_test(2);
		
	c(8)=-b_ub(3).* rx(4)+...
			 (1-p_test(4).*rx(4)-p_test(5).*rx(5)-...
				p_test(6).*rx(6)).*rx(3)./p_test(3);
			
	ceq=[];
elseif (type == 2) || (type == 3)
	% 这是K2和K3的（K2、K3可合并处理，因奖项数目一样）
	c(1)= 6e5- ...
			 (1-p_test(4).*rx(4)-p_test(5).*rx(5)-...
				p_test(6).*rx(6)-p_test(7).*rx(7)).*rx(1)./p_test(1);
	c(2)=-5e6+ ...
			 (1-p_test(4).*rx(4)-p_test(5).*rx(5)-...
				p_test(6).*rx(6)-p_test(7).*rx(7)).*rx(1)./p_test(1);
	c(3)= a_lb(1).* ...
			 (1-p_test(4).*rx(4)-p_test(5).*rx(5)-...
				p_test(6).*rx(6)-p_test(7).*rx(7)).*rx(2)./p_test(2)-...
			 (1-p_test(4).*rx(4)-p_test(5).*rx(5)-...
				p_test(6).*rx(6)-p_test(7).*rx(7)).*rx(1)./p_test(1);
	c(4)= a_lb(2).* ...
			 (1-p_test(4).*rx(4)-p_test(5).*rx(5)-...
				p_test(6).*rx(6)-p_test(7).*rx(7)).*rx(3)./p_test(3)-...
			 (1-p_test(4).*rx(4)-p_test(5).*rx(5)-...
				p_test(6).*rx(6)-p_test(7).*rx(7)).*rx(2)./p_test(2);
		
	c(5)= a_lb(3).* rx(4)-...
			 (1-p_test(4).*rx(4)-p_test(5).*rx(5)-...
				p_test(6).*rx(6)-p_test(7).*rx(7)).*rx(3)./p_test(3);

	c(6)=-b_ub(1).* ...
			 (1-p_test(4).*rx(4)-p_test(5).*rx(5)-...
				p_test(6).*rx(6)-p_test(7).*rx(7)).*rx(2)./p_test(2)+...
			 (1-p_test(4).*rx(4)-p_test(5).*rx(5)-...
				p_test(6).*rx(6)-p_test(7).*rx(7)).*rx(1)./p_test(1);
		
	c(7)=-b_ub(2).* ...
			 (1-p_test(4).*rx(4)-p_test(5).*rx(5)-...
				p_test(6).*rx(6)-p_test(7).*rx(7)).*rx(3)./p_test(3)+...
			 (1-p_test(4).*rx(4)-p_test(5).*rx(5)-...
				p_test(6).*rx(6)-p_test(7).*rx(7)).*rx(2)./p_test(2);
		
	c(8)=-b_ub(3).* rx(4)+...
			 (1-p_test(4).*rx(4)-p_test(5).*rx(5)-...
				p_test(6).*rx(6)-p_test(7).*rx(7)).*rx(3)./p_test(3);
	if (rx(7) == 0)
		c(9) = -1;
		c(10)= -1;
	else
		c(9) = a_lb(6) .* rx(7) - rx(6);
		c(10)=-b_ub(6) .* rx(7) + rx(6);
	end
	ceq=[];
elseif (type == 4)
	% K4型
	c(1)= 6e5- ...
			 (1-p_test(4).*rx(4)-p_test(5).*rx(5)).*rx(1)./p_test(1);
	c(2)=-5e6+ ...
			 (1-p_test(4).*rx(4)-p_test(5).*rx(5)).*rx(1)./p_test(1);
	c(3)= a_lb(1).* ...
			 (1-p_test(4).*rx(4)-p_test(5).*rx(5)).*rx(2)./p_test(2)-...
			 (1-p_test(4).*rx(4)-p_test(5).*rx(5)).*rx(1)./p_test(1);
	c(4)= a_lb(2).* ...
			 (1-p_test(4).*rx(4)-p_test(5).*rx(5)).*rx(3)./p_test(3)-...
			 (1-p_test(4).*rx(4)-p_test(5).*rx(5)).*rx(2)./p_test(2);
		
	c(5)= a_lb(3) .*rx(4)-...
			 (1-p_test(4).*rx(4)-p_test(5).*rx(5)).*rx(3)./p_test(3);

	c(6)=-b_ub(1).* ...
			 (1-p_test(4).*rx(4)-p_test(5).*rx(5)).*rx(2)./p_test(2)+...
			 (1-p_test(4).*rx(4)-p_test(5).*rx(5)).*rx(1)./p_test(1);
		
	c(7)=-b_ub(2).* ...
			 (1-p_test(4).*rx(4)-p_test(5).*rx(5)).*rx(3)./p_test(3)+...
			 (1-p_test(4).*rx(4)-p_test(5).*rx(5)).*rx(2)./p_test(2);
		
	c(8)=-b_ub(3) .* rx(4)+...
			 (1-p_test(4).*rx(4)-p_test(5).*rx(5)).*rx(3)./p_test(3);
	ceq=[];
else
	error('Error in function nonlcon!')
end

总结

  1. 决策论方法都需要用规划模型。彩票问题的目标比较难找，需要将心理效用和经济效用定量描述。
  2. 建模过程首先确定问题所属类型，本题就是最优决策问题。用目标规划方法比较合适。
  3. 确定目标(不便于量化，可利用心理、物理方法加以确定)
  4. 约束条件——明显约束以及隐含约束，隐含约束可以限制搜索范围。
  5. 模型求解——matlab(自主编程、基于判断)；lingo(整数规划)；mathematica(符号计算)。
  6. 确定求解工具后，就要确定求解算法，最原始的求解规划模型的方法是遍历算法，其次考虑智能算法提高求解效率。

2露天卡车调度问题2003B

规划类问题；模型容易建立；求解困难；求解策略：分布求解、逐级优化。

文档：2003B.note
链接：http://note.youdao.com/noteshare?id=23e7fbdb69be60ab5cbea2e70680e596&sub=D3E871044979486BBE00FA49F2BF5E71

c=[5.26	5.19	4.21	4	2.95	2.74	2.46	1.9	0.64	1.27	1.9	0.99	1.9	1.13	1.27	2.25	1.48	2.04	3.09	3.51	4.42	3.86	3.72	3.16	2.25	2.81	0.78	1.62	1.27	0.5	5.89	5.61	5.61	4.56	3.51	3.65	2.46	2.46	1.06	0.57	0.64	1.76	1.27	1.83	2.74	2.6	4.21	3.72	5.05	6.1];
A=[0.0154	0.0154	0.0154	0.0154	0.0154	0.0154	0.0154	0.0154	0.0154	0.0154	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
-0.0154	-0.0154	-0.0154	-0.0154	-0.0154	-0.0154	-0.0154	-0.0154	-0.0154	-0.0154	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0.0154	0.0154	0.0154	0.0154	0.0154	0.0154	0.0154	0.0154	0.0154	0.0154	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	-0.0154	-0.0154	-0.0154	-0.0154	-0.0154	-0.0154	-0.0154	-0.0154	-0.0154	-0.0154	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0.0154	0.0154	0.0154	0.0154	0.0154	0.0154	0.0154	0.0154	0.0154	0.0154	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	-0.0154	-0.0154	-0.0154	-0.0154	-0.0154	-0.0154	-0.0154	-0.0154	-0.0154	-0.0154	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
-0.005	-0.025	-0.015	0.015	0.005	0.025	0.015	0.005	0.025	0.005	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
-0.0104	0.0096	-0.0004	-0.0304	-0.0204	-0.0404	-0.0304	-0.0204	-0.0404	-0.0204	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	-0.005	-0.025	-0.015	0.015	0.005	0.025	0.015	0.005	0.025	0.005	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	-0.015	0.005	-0.005	-0.035	-0.025	-0.045	-0.035	-0.025	-0.045	-0.025	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	-0.005	-0.025	-0.015	0.015	0.005	0.025	0.015	0.005	0.025	0.005	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	-0.015	0.005	-0.005	-0.035	-0.025	-0.045	-0.035	-0.025	-0.045	-0.025	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
1	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	1	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	1	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	1	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	1	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	1	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	1	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	1	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	1	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	1
1	1	1	1	1	1	1	1	1	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	1	1	1	1	1	1	1	1	1	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1	1	1	1	1	1	1	1	1	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1	1	1	1	1	1	1	1	1	1	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1	1	1	1	1	1	1	1	1	1
1	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0
0	1	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0
0	0	1	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0
0	0	0	1	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0
0	0	0	0	1	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0
0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	1	0	0	0	0
0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	1	0	0	0
0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	1	0	0
0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	1	0
0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	1
];
b=[1.2154	-1.2	1.3154	-1.3	1.3154	-1.3	0	0	0	0	0	0	61.68831169	68.18181818	64.93506494	68.18181818	71.42857143	81.16883117	68.18181818	84.41558442	87.66233766	81.16883117	81.16883117	71.42857143	87.66233766	68.18181818	74.67532468	87.66233766	68.18181818	74.67532468	87.66233766	81.16883117	160	160	160	160	160	96	96	96	96	96	96	96	96	96	96]';
Aeq=[];
beq=[];
lb=zeros(1,50);
for i=1 :50
    ub(i)=inf;
end
[x,z]=linprog(c,A,b,Aeq,beq,lb,ub)

总结

  1. 思路：定目标——抽约束——求解。
  2. 求解技巧：循序渐进、分步求解、逐级优化，得到近似最优解。

3奥运会商圈规划问题2004A

应用了规划问题的求解方法、模拟退火算法和遗传算法。

文档：2004A.note
链接：http://note.youdao.com/noteshare?id=015712e03d5867ddd36b67fecf128e8f&sub=3191E37676B7428E891BB3FBB02E9057

%大商圈的布局问题——模拟退火算法
%MYFSAPLP SOLVE THE PLP IN THE THIRD PROBLEM BY SA ALGORITHM.
%首先只考虑n为偶数的情况
%变量初始化
clc
clf
clear;
xmin = 0; xmax = 200; ymin = 0; ymax = 150;
r = 25; n =6;
fval_every=1; fval_best=fval_every; fval_pro=fval_every;
lamdao=1e30; fs_every=1;
t0=98; tf=3; a=0.98; t=t0; t_mid=50;
p=1;
dfo=0;
while(t>tf)
	if p==1
		%产生新解
		for i=1:n./2
            x(i.*2-1)=r+2*r*(i-1);
            x(i*2)=r+2*r*(i-1);
            y(2*i-1)=r;
            y(2*i)=3*r;
		end
        for i=1:(n-1)
            for j=(i+1):n
                fs_every=fs_every.*sqrt((x(i)-x(j)).^2+(y(i)-y(j)).^2);
            end
        end
        fs_every=fs_every.^(1./n);
        fval_every=fs_every;
        dfval=fval_pro-fval_every;
		if dfval<0
            if p~=1
                x(i)=x0;
                y(i)=y0;
                fval_pro=fval_every;
                fs_pro=fs_every;
                dfo_pro=dfo;
            else
                fval_pro=fval_every;
                fs_pro=fs_every;
                dfo_pro=dfo;
            end
		else
            rand_jud=rand;
            if rand_jud>exp(-dfval./t)
                x(i)=x0;
                y(i)=y0;
                fval_pro=fval_every;
                fs_pro=fs_every;
                dfo_pro=dfo;
            else
                fs_every=fs_pro;
                fval_every=fval_pro;
                dfo=dfo_pro;
            end
		end
        if fval_every>fval_best
            x_best=x;
            y_best=y;
            fval_best=fval_every;
            fs_best=fs_every;
            dfo_bestk=dfo;
        end
        p=p-1;
        
	else
%         keyboard
        for k=1:500.*n
            i=ceil(rand.*n);
            x0=(xmin+r)+rand*((xmax-r)-(xmin+r));
            y0=(ymin+r)+rand*((ymax-r)-(ymin+r));
            for j=1:n
                if j~=i
                    dis0=sqrt((x(j)-x(i)).^2+(y(j)-y(i)).^2);
                    dis1=sqrt((x(j)-x0).^2+(y(j)-y0).^2);
                    if dis0 < 2*r
                        dfo=dfo-1;  
                    else
                        dfo=dfo+1;
                    end
                    if dis1 < 2*r
                        dfo=dfo+1;
                    else
                        dfo=dfo-1;
                    end
                    fs_every=fs_every./dis0.^(1./n);
                    fs_every=fs_every.*dis1.^(1./n);
                end
            end
            dfo_tmp=dfo./2;
            fval_every=fs_every-dfo_tmp.*lamdao./t;      
            
			dfval=fval_pro-fval_every;
			if dfval<0
                if p~=1
                    x(i)=x0;
                    y(i)=y0;
                    fval_pro=fval_every;
                    fs_pro=fs_every;
                    dfo_pro=dfo;
                else
                    fval_pro=fval_every;
                    fs_pro=fs_every;
                    dfo_pro=dfo;
                end
			else
                rand_jud=rand;
                if rand_jud>exp(-dfval./t)
                    x(i)=x0;
                    y(i)=y0;
                    fval_pro=fval_every;
                    fs_pro=fs_every;
                    dfo_pro=dfo;
                else
                    fs_every=fs_pro;
                    fval_every=fval_pro;
                    dfo=dfo_pro;
                end
			end
            if fval_every>fval_best
                x_best=x;
                y_best=y;
                fval_best=fval_every;
                fs_best=fs_every;
                dfo_bestk=dfo;
            end
        end
    end
    t=t.*a;
end

%求解结束
x_center = x_best;
y_center = y_best;
disp('圆心坐标分别为：')
circle_center = [x_center; y_center]';
disp(circle_center)
%绘图
hold on
for i = 1 : length(x_center)
    x_plot_tmp = linspace(x_center(i)-r,x_center(i)+r);
    y_plot_tmp_up = ...
        sqrt(r.^2 - (x_plot_tmp - x_center(i)).^2) ...
            + y_center(i);
    y_plot_tmp_down = ...
        - sqrt(r.^2 - (x_plot_tmp - x_center(i)).^2) ...
            + y_center(i);
    plot(x_plot_tmp,y_plot_tmp_up);
    plot(x_plot_tmp,y_plot_tmp_down);
end


%商圈的布局问题——遗传算法
% GA main
% using toolbox of GAOT version 5
clc
clear
bounds=ones(12,2);

global r

xmin=0; ymin=0; xmax=200; ymax=150;
r=25; n=6;
bounds(:,1)=zeros(12,1) + r;
bounds(1:6,2)=ones(6,1).*xmax - r;
bounds(7:12,2)=ones(6,1).*ymax - r;

% [x,endPop] = ga(bounds,'myfGAPLP',[],[],[1e-6 1 1]);
% [x,endPop] = ga('myfGAPLP',bounds,[],[],[1e-6 1 1]);
[x,endPop]= ga(@fGAPLP,12,[],[],[],[],bounds(:,1),bounds(:,2))
%subplot(1,2,1)
myfplotcircleGA(x,r,xmax,ymax);
xlabel('横坐标/m')
ylabel('纵坐标/m')
title('遗传算法的结果')



%fGAPLP.m
function fval=fGAPLP(x)
% options is required by the format of GAOT ver 5.

fval=1; dfo=0;
r=25;

n=6;lamdao=1e30;
for i=1:(n-1)
    for j=(i+1):n
        rtmp=sqrt((x(i)-x(j)).^2+(x(n+i)-x(n+j)).^2);
        if rtmp > 2*r
             fval=fval*rtmp;
             dfo=dfo+1;
        end
    end
end
fval=fval-dfo*lamdao;



%myfplotcircleGA.m
%% 以下为辅助函数的程序，放于独立的的m文件中
function y=myfplotcircleGA(xx,rtmp,xmax,ymax)
% plot circle

xtmp=xx(1:6);
ytmp=xx(7:12);
numlen=length(xtmp);
numsize=200;
t=linspace(0,2.*pi,numsize);
xplot=zeros(numsize,1);
yplot=zeros(numsize,1);
hold on
for j=1:numlen
	for i=1:numsize
		xplot(i)=xtmp(j)+rtmp*cos(t(i));
		yplot(i)=ytmp(j)+rtmp*sin(t(i));
    end
    plot(xplot,yplot)
 
end
y=[];
xlim([0 xmax]);
ylim([0 ymax]);

总结

  1. 综合的效益目标函数，包括经济效益、潜在效益、就业效益
  2. 计算机分析和实际情况结合给出了模型参数。通过lingo求解，得到了符合现实的结果。
  3. 运用商圈和雷利法则，对商区的MS布局给出了规划。
  4. 建模思路：定目标——抽约束——求解。
  5. 对于数据量的问题，处理、分析、统计数据，将数据变为有用的信息。
  6. 抽象模型目标时，考虑多些指标，因素，最后分析指标权重。
  7. 注意模型拓展：确定了各个商区超市数量以后，对超市的布局规划进行优化研究。

4交巡警服务平台的设置与调度2011B

处理大量数据、经典的最短路算法。

文档：2011B.note
链接：http://note.youdao.com/noteshare?id=722db129b099970671802c56bb1f44a6&sub=25357F9D3AAD45A7ADCEF83F19EC7311

题1：优化警台管辖范围

%求解直观上用0-1整数规划
%针对该问题对Floyd算法进行处理，即为各路口找出距其最近的静态


%链接：%https://pan.baidu.com/s/1lGmW12NXwzoEYCrSDQ12eQ 
%提取码：h3e2



%main.m
%% 最短路径的求解，以及现有的交巡警服务平台设置方案的合理性的判断
clc; clear all;
%% 数据准备
% 全市交通网络中路口节点节点坐标
A=xlsread('2011B_Table.xls', '全市交通路口节点数据','A2:C583'); 
% 全市交通路线
B=xlsread('2011B_Table.xls', '全市交通路口的路线','A2:B929');
% A区巡警台设置位置
PS_A=xlsread('2011B_Table.xls', '全市交巡警平台','B2:B21');

%% 计算最短距离矩阵
arn=size(A,1);
brn=size(B,1);
% 构建距离矩阵
D=ones(arn);
D(:)=inf;
for i=1:arn
    D(i,i)=0;
end

for i=1:brn
    m=B(i,1);   % 起点标号
    n=B(i,2);   % 终点标号
    d=sqrt((A(m,2)-A(n,2))^2+(A(m,3)-A(n,3))^2);
    D(m,n)=d;
    D(n,m)=d;
end
[dmin,path]=floyd(D);

%% 划分A区交巡警服务平台管辖范围
PN_A=A(1:92,1);
psn=size(PS_A, 1);
for i=1:92
    SR(i,1)=i;
    SR(i,2)=0;
    SR(i,3)=inf;
    SR(i,4)=0;
    for j=1:psn
        if SR(i,3)>dmin(i,j)
         SR(i,3)=dmin(i,j);
         SR(i,2)=j;
        end
    end
    if SR(i,3)*100>3000
       SR(i,4)=1;
    end    
end
%% 输出结果 
disp('A区交巡警服务平台管辖范围划分方案为：')
SR

     
     
% floyd.m
function [d,path]=floyd(a) 
% floyd   - 最短路问题  
%    a    - 距离矩阵是指i到j之间的距离可以是有向的   
%    d    - 最短路的距离  
%    path - 最短路的路径  
[n,lie]=size(a);   %  n为a的行数
d=a; 
path=zeros(n,n); 
% for i=1:n 
%    for j=1:n
%        if d(i,j)~=inf 
%            path(i,j)=j; %j是i的后续点 
%        end 
%    end 
% end  
for k=1:n
    for i=1:n 
      for j=1:n 
         if d(i,j)>d(i,k)+d(k,j) 
            d(i,j)=d(i,k)+d(k,j); 
            path(i,j)=path(i,k); 
         end 
      end 
   end 
end

题3：计算最佳新增警台数

%% 最短路径的求解，以及现有的交巡警服务平台设置方案的合理性的判断
clc; clear all; close all
%% 数据准备
% 全市交通网络中路口节点节点坐标
A=xlsread('2011B_Table.xls', '全市交通路口节点数据','A2:C583'); 
% 全市交通路线
B=xlsread('2011B_Table.xls', '全市交通路口的路线','A2:B929');
% A区巡警台设置位置
%PS_A=xlsread('2011B_Table.xls', '全市交巡警平台','B2:B21');

%% 计算最短距离矩阵
arn=size(A,1);
brn=size(B,1);
% 构建距离矩阵
D=ones(arn);
D(:)=inf;
for i=1:arn
    D(i,i)=0;
end

for i=1:brn
    m=B(i,1);   % 起点标号
    n=B(i,2);   % 终点标号
    d=sqrt((A(m,2)-A(n,2))^2+(A(m,3)-A(n,3))^2);
    D(m,n)=d;
    D(n,m)=d;
end
[dmin,path]=floyd(D);

%% 确定最佳新增警台数
alt_sta=21:92;
ppsl=1:20;
min_n=ones(1,5)*6;
% %% 启用双核加速计算(多核计算机才可以采用)
% if parpool('local', 6) == 0 
%     parpoolopen 6
% end
if isempty(gcp('nocreate'))
    parpool;
end
%% 计算各新增警台数与对应的不满足3分钟出警的路口数
for apsn=2:5
    com_pat=nchoosek(alt_sta,apsn);
    com_num=size(com_pat, 1);
    tmin_un_num=6;
    for k1=1:com_num
         npsl=[ppsl, com_pat(k1,:)]; % 更新警台设置方案
        %% 计算还有多少路口不满足3分钟出警 
         SR=zeros(92,4);
         for i=1:92
             i;
         SR(i,1)=i;
         SR(i,2)=0;
         SR(i,3)=inf;
         SR(i,4)=0;
                for j=1:(20+apsn)
                   if SR(i,3)>dmin(i,npsl(j))
                      SR(i,3)=dmin(i,npsl(j));
                      SR(i,2)=npsl(j);
                   end
                end
                if SR(i,3)*100>3000
                   SR(i,4)=1;
                end    
         end
         tun_num=sum(SR(:,4));
         if tun_num<tmin_un_num
            tmin_un_num=tun_num;
         end
         if tmin_un_num==0
             break
         end
         
    end
    min_n(1,apsn)=tmin_un_num;
end
     
% %% 关闭多核计算
% if parpool('local', 2) > 0%2个核心
%     parpool close
% end
% delete(gcp('nocreate'))
%% 显示结果
apn=[0, 2, 3, 4, 5]';
f_result=[apn, min_n']
figure(1) 
bar(apn,min_n)
xlabel('新增交警平台的数量');
ylabel('不满足3分钟出警的路口数');
title('不满足三分钟出警的路口数与新增服务台数的关系图','fontsize',12);
for i=1:4
    FT(1,i)=i+2*min_n(1,i+1);
end
opn=find(FT==min(FT));
disp(['最佳新增交警平台书:' num2str(opn+1)]);
T=2:5;
figure(2)
plot(T, FT,'-ko', 'LineWidth', 2);
xlabel('新增交警平台的数量');
ylabel('目标函数');
title('目标函数与新增交警平台的数关系图','fontsize',12);





% floyd算法的函数文件
function [d,path]=floyd(a) 
% floyd   - 最短路问题  
%    a    - 距离矩阵是指i到j之间的距离可以是有向的   
%    d    - 最短路的距离  
%    path - 最短路的路径  
[n,lie]=size(a);   %  n为a的行数
d=a; 
path=zeros(n,n); 
for i=1:n 
   for j=1:n
       if d(i,j)~=inf 
           path(i,j)=j; %j是i的后续点 
       end 
   end 
end  
for k=1:n
    for i=1:n 
      for j=1:n 
         if d(i,j)>d(i,k)+d(k,j) 
            d(i,j)=d(i,k)+d(k,j); 
            path(i,j)=path(i,k); 
         end 
      end 
   end 
end

题4：确定新增服务台位置(优化位置)

%main.m

%% 最短路径的求解，以及现有的交巡警服务平台设置方案的合理性的判断
clc; clear all; close all
%% 数据准备
% 全市交通网络中路口节点节点坐标
A=xlsread('2011B_Table.xls', '全市交通路口节点数据','A2:C583'); 
% 全市交通路线
B=xlsread('2011B_Table.xls', '全市交通路口的路线','A2:B929');
% 每个路口发案率(次数)
owl=xlsread('2011B_Table.xls', '全市交通路口节点数据','E2:E93');
%% 计算最短距离矩阵
arn=size(A,1)
brn=size(B,1);
% 构建距离矩阵
D=ones(arn);
D(:)=inf;
for i=1:arn
    D(i,i)=0;
end
for i=1:brn
    m=B(i,1);   % 起点标号
    n=B(i,2);   % 终点标号
    d=sqrt((A(m,2)-A(n,2))^2+(A(m,3)-A(n,3))^2);
    D(m,n)=d;
    D(n,m)=d;
end
[dmin,path]=floyd(D);

%% 确定最佳新增警台数
alt_sta=21:92;
ppsl=1:20;
min_n=ones(1,5)*6;
% %% 启用双核加速计算(多核计算机才可以采用)
% if matlabpool('size') == 0 
%     matlabpool open 2
% end
if isempty(gcp('nocreate'))
    parpool;
end

%% 计算各新增警台数与对应的不满足3分钟出警的路口数
    apsn=4;
    com_pat=nchoosek(alt_sta,apsn);
    com_num=size(com_pat, 1);
    avn=0;
    bfair_f=inf;
    for k1=1:com_num
         npsl=[ppsl, com_pat(k1,:)]; % 更新警台设置方案
        % 计算还有多少路口不满足3分钟出警 
         SR=zeros(92,4);
         for i=1:92
             i;
         SR(i,1)=i;
         SR(i,2)=0;
         SR(i,3)=inf;
         SR(i,4)=0;
                for j=1:(20+apsn)
                   if SR(i,3)>dmin(i,npsl(j))
                      SR(i,3)=dmin(i,npsl(j));
                      SR(i,2)=npsl(j);
                   end
                end
                if SR(i,3)*100>3000
                   SR(i,4)=1;
                end    
         end
         tun_num=sum(SR(:,4));
         
         if tun_num==0
          avn=avn+1;
          workload=zeros(1, 24);
               for j=1:24
                 for u=1:92
                     if SR(u,2)==npsl(1,j)
                     workload(1,j)=workload(1,j)+owl(u,1);
                     end
                  end
               end
               fair_f=var(workload,1);
               fair_r(avn,:)=[avn, fair_f,  com_pat(k1,:)];             
               
              if fair_f<bfair_f
                   bfair_f=fair_f;
                   ba=com_pat(k1,:);
                   bs=SR;
              end
         end       
    end      
% %% 关闭多核计算
% if matlabpool('size') > 0
%     matlabpool close
% end
%% 显示结果
abpn=0;
for i=1:avn
    if fair_r(i,2)==bfair_f
        abpn=abpn+1;
        abps(abpn,:)=[abpn, fair_r(i,2:end)];
    end
end
disp(['所有可行方案数' num2str(avn)]);
disp(['最佳方案数' num2str(abpn)]);
disp('所有最佳方案:')
abps
disp('最佳方案对应的工作均衡度:')
bfair_f
plot(fair_r(:,1),fair_r(:,2),'-ko', 'LineWidth', 2);
set(gca,'linewidth',2) ;
xlabel('可行方案编号'); ylabel('工作均衡度');
title('各可行方案的工作均衡度','fontsize',12);       



%floyd.m
% floyd算法的函数文件
function [d,path]=floyd(a) 
% floyd   - 最短路问题  
%    a    - 距离矩阵是指i到j之间的距离可以是有向的   
%    d    - 最短路的距离  
%    path - 最短路的路径  
[n,lie]=size(a);   %  n为a的行数
d=a; 
path=zeros(n,n); 
for i=1:n 
   for j=1:n
       if d(i,j)~=inf 
           path(i,j)=j; %j是i的后续点 
       end 
   end 
end  
for k=1:n
    for i=1:n 
      for j=1:n 
         if d(i,j)>d(i,k)+d(k,j) 
            d(i,j)=d(i,k)+d(k,j); 
            path(i,j)=path(i,k); 
         end 
      end 
   end 
end

总结

  1. 目标规划问题
  2. 确定决策变量——甄选目标函数——抽象约束条件
  3. 模型求解：最短路算法Floyd算法

5葡萄酒评价2012A

文档：2012A.note
链接：http://note.youdao.com/noteshare?id=7410452f1b1b54a5812c1bf54b49ffb1&sub=256456C2802D41E2AC2E9B2C71A9E543

题1：判断评价结果有无显著性差异

%题1,2链接：https://pan.baidu.com/s/1JEJOHzXwM3Uw3N0uWnm9iw 
%提取码：j4na
%题4链接：https://pan.baidu.com/s/1c5cKZh_WplpZpKKwSJXejA 
%提取码：yixo


%酒样得分T检查
%% 2012A_question1_T evaluation
% 程序15-1
%--------------------------------------------------------------------------
%% 数据准备
% 清空环境变量
clear all
clc

%导入数据
X1=xlsread('2012A_T1_processed.xls', 'T1_red_grape', 'D3:M272');  
X2=xlsread('2012A_T1_processed.xls', 'T2_red_grape', 'D3:M272');
X3=xlsread('2012A_T1_processed.xls', 'T1_white_grape', 'D3:M282');
X4=xlsread('2012A_T1_processed.xls', 'T2_white_grape', 'D3:M282');

%% 红葡萄酒T检验计算过程
[m1,n1]=size(X1);
K1=27; 
% 计算每个样品的总得分
for i=1:K1
      for j=1:n1
      SX1(i,j)=sum(X1(10*i-9:10*i,j));
      SX2(i,j)=sum(X2(10*i-9:10*i,j));
      end
end
% 计算每组样品得分的均值
for i=1:K1
    Mean1(i)=mean(SX1(i,:));
    Mean2(i)=mean(SX2(i,:));
end
% 计算检验值
for i=1:K1
    S1(1,i)=(sum((SX1(i,:)-Mean1(i)).^2)+sum((SX2(i,:)-Mean2(i)).^2))/(n1*(n1-1));
    T1(1,i)=(Mean1(i)-Mean2(i))/(sqrt(S1(1,i)));
end
AT_R=abs(T1);  
M_AT_R=mean(AT_R);
%% 白葡萄酒T检验计算过程
[m2,n2]=size(X3);
K2=28;
% 计算每个样品的总得分
for i=1:K2
      for j=1:n2
      SX3(i,j)=sum(X3(10*i-9:10*i,j));
      SX4(i,j)=sum(X4(10*i-9:10*i,j));
      end
end
% 计算每组样品得分的均值
for i=1:K2
    Mean3(i)=mean(SX3(i,:));
    Mean4(i)=mean(SX4(i,:));
end
% 计算检验值
for i=1:K2
    S2(1,i)=(sum((SX3(i,:)-Mean3(i)).^2)+sum((SX4(i,:)-Mean4(i)).^2))/(n2*(n2-1));
    T2(1,i)=(Mean3(i)-Mean4(i))/(sqrt(S2(1,i)));
end
AT_W=abs(T2);  
M_AT_W=mean(AT_W);
%% 结果显示与比较
a=2.102; % T(0.05,2,18)=2.101
b=2.878; % T(0.01,2,18)=2.878
set(gca,'linewidth',2) 
% 红酒结果
for i=1:K1
    Ta1(i)=a;
    Tb1(i)=b;
end
t1=1:K1;
subplot(2,1,1);
plot(t1,AT_R,'*k-',t1,Ta1,'r-',t1,Tb1,'-.b', 'LineWidth', 2)
title('红酒显著性检验结果','fontsize',14)
legend('T检验值', 'T(0.05)值', 'T(0.01)值')
xlabel('样品号'), ylabel('T检验值')

% 白酒结果
for i=1:K2
    Ta2(i)=a;
    Tb2(i)=b;
end
t2=1:K2;
subplot(2,1,2);
plot(t2,AT_W,'*k-',t2,Ta2,'r-',t2,Tb2,'-.b', 'LineWidth', 2)
title('白酒显著性检验结果','fontsize',14)
legend('T检验值', 'T(0.05)值', 'T(0.01)值')
xlabel('样品号'), ylabel('T检验值')
% 显示平均检验结果
disp(['两组品酒师对红酒的平均显著性T检验值:' num2str(M_AT_R)]);
disp(['两组品酒师对白酒的平均显著性T检验值:' num2str(M_AT_W)]);

题2：评价结果稳定性

%两组品酒师评价
%% 2012A_question1_T evaluation
% 程序16-2
%--------------------------------------------------------------------------
%% 数据准备
% 清空环境变量
clear all
clc

%导入数据(白葡萄酒)
% X1=xlsread('2012A_T1_processed.xls', 'T1_white_grape', 'D3:M282');
% X2=xlsread('2012A_T1_processed.xls', 'T2_white_grape', 'D3:M282');
%导入数据(红葡萄酒)
X1=xlsread('2012A_T1_processed.xls', 'T1_red_grape', 'D3:M272');
X2=xlsread('2012A_T1_processed.xls', 'T2_red_grape', 'D3:M272');

%% 计算每组品酒师对每个样品的方差
[m,n]=size(X1);
K=27;
% K=28 白葡萄酒
% 计算每个样品的总得分
for i=1:K
      for j=1:n
      SX1(i,j)=sum(X1(10*i-9:10*i,j));
      SX2(i,j)=sum(X2(10*i-9:10*i,j));
      end
      u0(i)=mean([SX1(i,:), SX2(i,:)]);
end
% 计算方差
for i=1:K
    SD1(i,:)=(SX1(i,:)-u0(i)).*(SX1(i,:)-u0(i));
    SD2(i,:)=(SX2(i,:)-u0(i)).*(SX2(i,:)-u0(i));
end

%% 结果显示与比较
for i=1:K
    TSD(1,i)=sum(SD1(i,:));
    TSD(2,i)=sum(SD2(i,:));
end
t=1:K;
plot(t,TSD(1,:),'*k-',t,TSD(2,:),'ok--', 'LineWidth', 2)
set(gca,'linewidth',2);
legend('一组方差','二组方差')
xlabel('红葡萄酒样品编号'); ylabel('红葡萄酒评价方差');
TSD1=sum(TSD(1,:));
TSD2=sum(TSD(2,:));
disp(['一组对白葡萄酒总方差:' num2str(TSD1)]);
disp(['二组对白葡萄酒总方差:' num2str(TSD2)]);

题3：根据葡萄酒质量对葡萄分级(聚类)

%% 用聚类法确定葡萄酒分级
clc, clear all, close all
%% 需要聚类的数据
% 红葡萄酒质量评分数据
A=[79.95	75	80.45	78.15	76.25	71.95	75.85	71.85	76.65	77.05...
    71.85	67.84	69.9    74.55	75.4	70.65	79.55	74.9	74.3	77.2...
    77.8	75.2	76.65	74.7	78.3    77.8	70.9	80.45]; 
% 白葡萄酒质量评分数据
% A=[79.95	75	80.45	78.15	76.25	71.95	75.85	71.85	76.65	77.05...
% 71.85	67.85	69.9	74.55	75.4	70.65	79.55	74.9	74.3	77.2...
%     77.8	75.2	76.65	74.7	78.3	77.8	70.9]
%% 用k-Means聚类法确定最佳的聚类数
X=A';
numC=7;
for i=1:numC
    kidx = kmeans(X,i);
    silh = silhouette(X,kidx); %计算轮廓值
    silh_m(i) = mean(silh);    %计算平均轮廓值
end

figure
plot(1:numC,silh_m,'ko-', 'linewidth',2)
set(gca,'linewidth',2);
xlabel('类别数')
ylabel('平均轮廓值')
title(' 不同类别对应的平均轮廓值')

% 绘制2至5类时的轮廓值分布图
 figure
 set(gca,'linewidth',2);
for i=2:5
    kidx = kmeans(X,i);
    subplot(2,2,i-1);
    [~,h] = silhouette(X,kidx);
    set(gca,'linewidth',2);
    title([num2str(i), '类时的轮廓值 ' ])
    snapnow
    xlabel('轮廓值');
    ylabel('类别数');
end

%% K-means聚类过程，并将结果显示出来
[idx,ctr]=kmeans(A',4); % 用K-means法聚类
% 提取同一类别的样品号
c1=find(idx==1); c2=find(idx==2); 
c3=find(idx==3); c4=find(idx==4); 
figure
F1 = plot(find(idx==1), A(idx==1),'k:*', ...
     find(idx==2), A(idx==2),'k:o', ...
     find(idx==3), A(idx==3),'k:p', ...
     find(idx==4), A(idx==4),'k:d');
set(gca,'linewidth',2);
set(F1,'linewidth',2, 'MarkerSize',8);
xlabel('编号','fontsize',12);
ylabel('得分','fontsize',12);
title('Kmeans方法聚类结果')
disp('聚类结果：');
disp(['第1类:' ,'中心点：',num2str(ctr(1)),'  ','该类样品编号：', num2str(c1')]);
disp(['第2类:' ,'中心点：',num2str(ctr(2)),'  ','该类样品编号：', num2str(c2')]);
disp(['第3类:' ,'中心点：',num2str(ctr(3)),'  ','该类样品编号：', num2str(c3')]);
disp(['第4类:' ,'中心点：',num2str(ctr(4)),'  ','该类样品编号：', num2str(c4')]);

%% 层次聚类
X=A';
Y=pdist(X);  %计算样品间的欧式距离
Z=linkage(Y, 'average');  % 利用类平均法创建系统聚类树
cn=size(X);
clabel=1:cn;
clabel=clabel';
figure
F2 = dendrogram(Z); % 绘制聚类树形图
set(gca,'linewidth',2);
title('层次聚类法聚类结果')
set(F2,'linewidth',2);
ylabel('标准距离');
%% Fuzzy C-means聚类
X=A';
[center,U] = fcm(X,4);
Cid1 = find(U(1,:) ==max(U));
Cid2 = find(U(2,:) ==max(U));
Cid3 = find(U(3,:) ==max(U));
Cid4 = find(U(4,:) ==max(U));
figure
F3=plot(Cid1, A(Cid1),'k:*', ...
     Cid2, A(Cid2),'k:o', ...
     Cid3, A(Cid3),'k:p', ...
     Cid4, A(Cid4),'k:d');

set(gca,'linewidth',2);
set(F3,'linewidth',2, 'MarkerSize',8);
xlabel('编号');
ylabel('得分');
title('Fuzzy C-means方法聚类结果')
%% 2012年全国赛A题第二问求解示例程序

题4：根据葡萄酒理化指标分级(聚类)

为保持指标级别一致性，统一选择一级指标作为研究对象。

%% PCA-kmeans方法实现对葡萄的理化指标进行聚类
%--------------------------------------------------------------------------
%% 数据导入及处理
clc, clear all, close all
%  A=xlsread('2012A_Table2.xls','葡萄酒指标汇总', 'C3:J29');% 红葡萄酒
A=xlsread('2012A_Table2.xls','葡萄酒指标汇总', 'C33:J60');% 白葡萄酒

%  数据标准化处理
a=size(A,1);  
b=size(A,2);  
for i=1:b
    SA(:,i)=(A(:,i)-mean(A(:,i)))/std(A(:,i)); 
end

%% 计算相关系数矩阵的特征值和特征向量
CM=corrcoef(SA);  % 计算相关系数矩阵(correlation matrix)
[V, D]=eig(CM);  % 计算特征值和特征向量

for j=1:b
    DS(j,1)=D(b+1-j, b+1-j); % 对特征值按降序进行排序
end
for i=1:b
    DS(i,2)=DS(i,1)/sum(DS(:,1)); %贡献率
    DS(i,3)=sum(DS(1:i,1))/sum(DS(:,1)); %累积贡献率
end

%% 选择主成分及对应的特征向量
T=0.8;  % 主成分信息保留率.
for K=1:b
    if DS(K,3)>=T
        Com_num=K;
        break;
    end
end

% 提取主成分对应的特征向量
for j=1:Com_num
    PV(:,j)=V(:,b+1-j);
end

%%  计算各评价对象的主成分得分
new_score=SA*PV;
for i=1:a
    total_score(i,1)=sum(new_score(i,:));
    total_score(i,2)=i;
end
result_report=[new_score, total_score]; % 将各主成分得分与总分放在同一个矩阵中
result_report=sortrows(result_report,(K+2)); % 按总分降序排序

%% Kmeans聚类及结果报告
A=result_report(:,(K+1));
[idx,ctr]=kmeans(A,4);
[m,n]=size(A);
t1=ones(1,n)*30;
c1=find(idx==1); c2=find(idx==2); c3=find(idx==3); c4=find(idx==4);
h=plot(t1,A,'ko',c1,A(idx==1),'k--*', c2,A(idx==2),'k--s', c3,A(idx==3),'k--d', c4,A(idx==4),'k--p');
xlabel('白葡萄酒样品编号','fontsize',12);
ylabel('主成分得分','fontsize',12);
title('白葡萄酒理化指标聚类图','fontsize',12)
set(h, 'MarkerSize',8, 'MarkerFaceColor','k');
set(gca,'linewidth',2) ;

disp('主成分得分(最后1列为样本编号，倒数第2列为总分，前面为各主成分得分)')
result_report  
disp('分类结果：');
disp(['第1类:' ,'中心点：',num2str(ctr(1)),'  ','该类样品编号：', num2str(c1')]);
disp(['第2类:' ,'中心点：',num2str(ctr(2)),'  ','该类样品编号：', num2str(c2')]);
disp(['第3类:' ,'中心点：',num2str(ctr(3)),'  ','该类样品编号：', num2str(c3')]);
disp(['第4类:' ,'中心点：',num2str(ctr(4)),'  ','该类样品编号：', num2str(c4')]);

总结

  1. t检验、K-means、PCA

6经验

  1. 知识：图论、排队论、泛函、运筹学
  2. 书籍：姜启源、谢金星、叶其孝、朱道元
  3. 看书：不求甚解，要知道基本思路。
  4. 算法：遍历、二分、数据结构；智能优化算法(用于求解)
  5. set命令重设图形的句柄属性，可绘制精美的图形
  6. 好的数学模型：有效、简单、恰当地解决实际问题(易于理解)
  7. 好的求解方法：快速、准确给出最优解
  8. 高质量论文：20页论文
  9. 进度：1小时确定选题、第一天上午建好数学模型并确定求解方法。
  10. 选题后——各自看题——列出疑点——讨论——理解统一(40min)
  11. 建模思路——讨论——写第一问——matlab、mathematics求解。(下午两点之前)——第二问
  12. 第一天晚上10点，完成了所有建模和求解。
  13. 接下来的时间对算法改进，原来的问题上加入新课题。
  14. 目标建模法——如果没有思路，想象最后结果是什么，再分析实现目标的途径。
HTUer的编程之路
关注
  • 4
    点赞
  • 踩
  • 20
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 0
    评论
专栏目录
一个数学建模的目标规划问题：奥运会商圈规划问题
搬家啦x
02-05 5733
奥运会商圈规划问题 1 问题描述        2008奥运会期间，在比赛主场馆周边地区需要建设迷你超市（MS）网，以满足观众的购物需求。        基本要求：满足奥运会期间的购物需求、分布基本均衡和商业上盈利。   2 基本假设、名词约定及符号说明 2.1 基本假设（包含对目标要求的基本理解） （1）  各场馆相互独立； （2）  MS分布均衡是指各商区的MS的数量相等或近似
matlab开发-组合Nchoosek
08-25
matlab开发-组合Nchoosek。此函数查找n choose k的组合。
《Matlab在数学建模中的应用》笔记系列1-线性规划
Eajack's CSDN blog
09-02 433
###第一次在CSDN上写博客，暂时用于笔记用途###
MATLAB中nchoosek的用法
weixin_45763816的博客
10-26 6480
MATLAB中nchoosek的用法为： 得到二项式系数或所有的组合，读起来有点拗口。其实通俗一点就是从多个元素中选出一定数量的所有组合（不分先后）。 用法为一：a=nchoosek(n,k)，b,k均为非负整数，则a为二项式系数。 tips：二项式系数的意义是从n件物件中，不分先后地选取k件的方法总数，因此也叫做组合数。 用法二：b=nchoosek（v,k），v为向量，k为非负整数，则b为从向量v中一次取 k 个元素的所有组合。矩阵b 有 k 列和 m!/((m–k)! k!) 行，其中 m 为 向量v
nchoosek--二项式系数或所有的组合数
qq_18343569的博客
12-13 1718
【功能简介】求组合数。 【语法格式】 1．C=nchoosek(n,k) 参数n、k为非负整数，返回 =n!/((n-k)!k!)。 2．C=nchoosek(v,k) v是长度为n的行向量，返回的C是一个矩阵，矩阵的行包含从n个元素中取k个的所有可能组合。因此，C是nchoosek(n,k)行k列矩阵。 【实例5.30】求从1~5中取出3个数的所有不同取法。
Maltab在数学建模中的应用（第二版） 含源程序
08-29
Maltab在数学建模中的应用（第二版） 含源程序,书本程序实例详细，完整。很值得参加数模比赛的同学学习。
Maltab在数学建模中的应用（第二版）_卓金武等
10-26
《Maltab在数学建模中的应用（第二版）》是由卓金武等专家编著的一本关于如何利用MATLAB进行数学建模的专著。这本书深入浅出地介绍了MATLAB这一强大的数值计算和可视化工具在数学建模领域的广泛应用，为读者提供了...
Applicationmathematicalmodeling_Maltab在数学建模中的应用_gettingq2l_
10-03
《Matlab在数学建模中的应用——以"Getting q2l"为例》 数学建模是一种将实际问题抽象成数学模型的过程，它在科学研究、工程计算、经济预测等领域发挥着重要作用。而Matlab，作为一款强大的数值计算和可视化工具，...
Maltab在数学建模中的应用（第二版） 源程序.rar
08-20
《Matlab在数学建模中的应用（第二版）》是一本深入探讨如何使用Matlab进行数学建模的专著。该书的源程序文件集合，包含了一系列精心设计的实例，旨在帮助读者掌握Matlab在解决实际问题时的强大功能。源程序经过测试...
源程序_Maltab在数学建模中的应用（第二版）_卓金武等
09-30
《源程序_Maltab在数学建模中的应用（第二版）_卓金武等》是一部深入探讨Matlab在数学建模领域应用的专著。本书由卓金武等专家精心编写，旨在帮助读者掌握如何利用Matlab的强大功能解决实际的数学建模问题。通过实例...
nchoosek_enum:带有枚举选择的nchoosek。 此函数永远不会在内存中构建完整组合数组-matlab开发
05-28
此功能可用于处理基数k集的基数k的所有子集，而无需将其存储在内存中。 典型的电话如下 n = 10； k = 3； chunksize = 10; i = 1：块大小； norepetition =假； 虽然真实如果没有重复c = nchoosek_enum（n，k，i）; 别的c = nchoosek_enum（n + k-1，k，i）; c = c-（0：k-1）; 结尾对于r = 1：size（c，1） cr = c（r，:); endcomb = any（isnan（cr））; 如果端梳休息结尾％用cr做某事fprintf（'％s \ n'，mat2str（cr）） 结尾如果端梳休息结尾我=我+块大小; 结尾
nchoosek排列组合
06-03
该代码课用于排列组合问题的求解即从n个元素中取m个，m<=n的问题
matlab 遗传算法组合优化实例
08-16
matlab多目标遗传算法优化实例，注释详尽。序安装遗传算法工具箱。
基于K均值聚类的葡萄酒品种判别
学习ing的香干
06-21 8869
摘要 本文基于K均值聚类方法，首先利用使得组间平方和较大、组内平方和较小的方法确定聚类簇K，再对葡萄酒数据进行聚类，得到误判率为3.37%。为描述不同类别葡萄酒的特点，利用主成分分析，探究影响分类的因素，发现类别1与非类别1的差异主要是脯氨酸含量的高低。 关键词：K均值聚类；葡萄酒；主成分分析 ABSTRACT This article trys to classify distinguishe...
一个matlab里的排列组合或者说是循环嵌套的层数不定的问题
hank的主页
07-01 3614
我要完成类似这么一个工作：我有一个n维的坐标，[a,b,c....n]坐标的范围已知，如何将所有的坐标都填到一个矩阵里去，举例说明，假设有[a,b,c]取值在[3,3,3]到[4,4,4]之间，那么应该这之间的所有坐标为444,344,434,443,334,343,433,333也就是2*2*2个，这是一个排列组合的问题，可是我试了用perms好像无法实现我的要求，而如果用循环嵌套的方式，因为n...
matlab列出所有排列组合数
热门推荐
ego132639的博客
05-30 1万+
设有三类机组类型，台数分别为4，3，2台，每种机组的运行台数情况如下： x1=[0 1 2 3 4]; %表示类型一机组运行的台数可能为0台或1台或2台或3台或4台 x2=[0 1 2 3]; %表示类型一机组运行的台数可能为0台或1台或2台或3台 x3=[0 1 2]; %表示类型一机组运行的台数可能为0台或1台或2台 [x3,x2,x1] = ndgrid(x3,x2,x1); combi...
maltab在数学建模中的应用(第二版)_卓金武等
最新发布
01-20
《maltab在数学建模中的应用(第二版)》一书由卓金武等人撰写，详细介绍了maltab在数学建模中的广泛应用。maltab是一种强大的数学计算软件，被广泛用于工程、科学和数学领域，能够帮助研究人员快速、精确地进行数据...

“相关推荐”对你有帮助么？

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

HTUer的编程之路

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付：¥1
获取中
扫码支付

您的余额不足，请更换扫码支付或充值

打赏作者

实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值