数学建模 -- 蒙特卡洛算法

来自清风的数学建模课程，主要是用于自己复习看，所以截图较多

蒙特卡洛算法

定义

原理

大数定理可知，当样本容量足够大时，事件发生的频率就是概率

讨论

蒙特卡罗方法的示例

三门问题

%% （2）代码部分（在成功的条件下的概率）
n = 100000;  % n代表蒙特卡罗模拟重复次数
a = 0;  % a表示不改变主意时能赢得汽车的次数
b = 0;  % b表示改变主意时能赢得汽车的次数
for i= 1 : n  % 开始模拟n次
    x = randi([1,3]);  % 随机生成一个1-3之间的整数x表示汽车出现在第x扇门后
    y = randi([1,3]);  % 随机生成一个1-3之间的整数y表示自己选的门
    % 下面分为两种情况讨论：x=y和x~=y
    if x == y   % 如果x和y相同，那么我们只有不改变主意时才能赢
        a = a + 1;     b = b + 0;
    else  % x ~= y ，如果x和y不同，那么我们只有改变主意时才能赢
        a = a + 0;     b = b +1;
    end
end
disp(['蒙特卡罗方法得到的不改变主意时的获奖概率为：', num2str(a/n)]);
disp(['蒙特卡罗方法得到的改变主意时的获奖概率为：', num2str(b/n)]);


%% （3）考虑失败情况的代码(无条件概率)
n = 100000;  % n代表蒙特卡罗模拟重复次数
a = 0;  % a表示不改变主意时能赢得汽车的次数
b = 0;  % b表示改变主意时能赢得汽车的次数
c = 0;  % c表示没有获奖的次数
for i= 1 : n  % 开始模拟n次
    x = randi([1,3]);  % 随机生成一个1-3之间的整数x表示汽车出现在第x扇门后
    y = randi([1,3]);  % 随机生成一个1-3之间的整数y表示自己选的门
    change = randi([0, 1]); % change =0  不改变主意，change = 1 改变主意
    % 下面分为两种情况讨论：x=y和x~=y
    if x == y   % 如果x和y相同，那么我们只有不改变主意时才能赢
        if change == 0  % 不改变主意
        	a = a + 1; 
        else  % 改变了主意
            c= c+1;
        end
    else  % x ~= y ，如果x和y不同，那么我们只有改变主意时才能赢
         if change == 0  % 不改变主意
        	c = c + 1; 
        else  % 改变了主意
            b= b + 1;
         end
    end
end
disp(['蒙特卡罗方法得到的不改变主意时的获奖概率为：', num2str(a/n)]);
disp(['蒙特卡罗方法得到的改变主意时的获奖概率为：', num2str(b/n)]);
disp(['蒙特卡罗方法得到的没有获奖的概率为：', num2str(c/n)]);

模拟排队问题

%% （3）问题1的代码
clear
tic  %计算tic和toc中间部分的代码的运行时间
i = 1;  % i表示第i个客户，最开始取i=1
w = 0;  % w用来表示所有客户等待的总时间，初始化为0
e0 = 0;  c0 = 0;   % 初始化e0和c0为0
x(1) = exprnd(10);  % 第0个客户(假想的)和第1个客户到达的时间间隔
c(1) = c0 + x(1);  % 第1个客户到达的时间
b(1) = c(1); % 第1个客户的开始服务的时间
while b(i) <= 480  % 开始设置循环，只要第i个顾客开始服务的时间(时刻)小于480，就可以对其服务（银行每天工作8小时，折换为分钟就是480分钟）
    y(i) = normrnd(10,2); % 第i个客户的服务持续时间，服从均值为10方差为4(标准差为2)的正态分布
    if y(i) < 1  % 根据题目的意思：若服务持续时间不足一分钟，则按照一分钟计算
        y(i) = 1;
    end
    e(i) = b(i) + y(i); % 第i个客户结束服务的时间 = 第i个客户开始服务的时间 + 第i个客户的服务持续时间
    wait(i) = b(i) - c(i); % 第i个客户等待的时间 = 第i个客户开始服务的时间 - 第i个客户到达银行的时间
    w = w + wait(i); % 更新所有客户等待的总时间
    i = i + 1; % 增加一名新的客户
    x(i) = exprnd(10); % 这位新客户和上一个客户到达的时间间隔
    c(i) = c(i-1) + x(i); % 这位新客户到达银行的时间 = 上一个客户到达银行的时间 + 这位新客户和上一个客户到达的时间间隔
    b(i) = max(c(i),e(i-1)); % 这个新客户开始服务的时间取决于其到达时间和上一个客户结束服务的时间
end
n = i-1; % n表示银行一天8小时一共服务的客户人数
t = w/n; % 客户的平均等待时间
disp(['银行一天8小时一共服务的客户人数为: ',num2str(n)])
disp(['客户的平均等待时间为: ',num2str(t)])
toc  %计算tic和toc中间部分的代码的运行时间

%% （4）问题2的代码
clear
tic  %计算tic和toc中间部分的代码的运行时间
day = 100;  % 假设模拟100天
n = zeros(day,1); % 初始化用来保存每日接待客户数结果的矩阵
t = zeros(day,1); % 初始化用来保存每日客户平均等待时长的矩阵
for k = 1:day
    i = 1;  % i表示第i个客户，最开始取i=1
    w = 0;  % w用来表示所有客户等待的总时间，初始化为0
    e0 = 0;  c0 = 0;   % 初始化e0和c0为0
    x(1) = exprnd(10);  % 第0个客户(假想的)和第1个客户到达的时间间隔
    c(1) = c0 + x(1);  % 第1个客户到达的时间
    b(1) = c(1); % 第1个客户的开始服务的时间
    while b(i) <= 480  % 开始设置循环，只要第i个顾客开始服务的时间(时刻)小于480，就可以对其服务（银行每天工作8小时，折换为分钟就是480分钟）
        y(i) = normrnd(10,2); % 第i个客户的服务持续时间，服从均值为10方差为4(标准差为2)的正态分布
        if y(i) < 1  % 根据题目的意思：若服务持续时间不足一分钟，则按照一分钟计算
            y(i) = 1;
        end
        e(i) = b(i) + y(i); % 第i个客户结束服务的时间 = 第i个客户开始服务的时间 + 第i个客户的服务持续时间
        wait(i) = b(i) - c(i); % 第i个客户等待的时间 = 第i个客户开始服务的时间 - 第i个客户到达银行的时间
        w = w + wait(i); % 更新所有客户等待的总时间
        i = i + 1; % 增加一名新的客户
        x(i) = exprnd(10); % 这位新客户和上一个客户到达的时间间隔
        c(i) = c(i-1) + x(i); % 这位新客户到达银行的时间 = 上一个客户到达银行的时间 + 这位新客户和上一个客户到达的时间间隔
        b(i) = max(c(i),e(i-1)); % 这个新客户开始服务的时间取决于其到达时间和上一个客户结束服务的时间
    end
    n(k) = i-1; % n(k)表示银行第k天服务的客户人数
    t(k) = w/n(k); % t(k)表示该银行第k天客户的平均等待时间
end
disp([num2str(day),'个工作日中，银行每日平均服务的客户人数为: ',num2str(mean(n))])
disp([num2str(day),'个工作日中，银行每日客户的平均等待时间为: ',num2str(mean(t))])
toc  %计算tic和toc中间部分的代码的运行时间

有约束的非线性规划问题

求得是近似解 局部最优解

0-1规划

旅行商问题