title: 蒙特卡羅（洛）模擬
date: 2020-02-27 21:26:53
categories: 數學建模
tags: [ MATLAB, 模擬]
mathjax: true
引例
布豐投針實驗
法國數學家布豐（1707-1788）最早設計了投針試驗。
這一方法的步驟是：
1） 取一張白紙，在上面畫上許多條間距為a的平行線。
2） 取一根長度為l（l≤a） 的針，隨機地向畫有平行直線的紙上擲n次，觀察針與直線相交的次數，記為m。
3）計算針與直線相交的概率．
18世紀，法國數學家布豐提出的“投針問題”，記載於布豐1777年出版的著作中：“在平面上畫有一組間距為a的平行線，將一根長度為l（l≤a）的針任意擲在這個平面上，求此針與平行線中任一條相交的概率。”
布豐本人證明了，這個概率是：
\[p=\frac{2l}{\pi a} \]
(其中π為圓周率）
由於它與π有關，於是人們想到利用投針試驗來估計圓周率的值。
布豐驚奇地發現：有利的扔出與不利的扔出兩者次數的比，是一個包含π的表示式．如果針的長度等於a/2，那么有利扔出的概率為1/π．扔的次數越多，由此能求出越為精確的π的值。
下面是利用這個公式，用概率的方法得到圓周率的近似值的一些資料。
| 試驗者 | 時間 | 投擲次數 | 相交次數 | 圓周率估計值 |
|---|---|---|---|---|
| Wolf | 1850年 | 5000 | 2532 | 3.1596 |
| Smith | 1855年 | 3204 | 1218.5 | 3.1554 |
| C.De Morgan | 1860年 | 600 | 382.5 | 3.137 |
| Fox | 1884年 | 1030 | 489 | 3.1595 |
| Lazzerini | 1901年 | 3408 | 1808 | 3.1415929 |
| Reina | 1925年 | 2520 | 859 | 3.1795 |
證明
代碼
%% 蒙特卡羅用於布豐投針實驗
%% （1）預備知識
% rand(m,n)函數產生由在[0,1]之間均勻分布的隨機數組成的m行n列的矩陣（或稱為數組）。
rand(5,4)
% 0.8300 0.1048 0.2396 0.4398
% 0.5663 0.1196 0.8559 0.5817
% 0.9281 0.2574 0.3013 0.9355
% 0.3910 0.3173 0.2108 0.1676
% 0.3645 0.4372 0.8819 0.9232
rand(3) % 若只給一個輸入，則會生成一個方陣
% 0.1709 0.4951 0.0541
% 0.9374 0.8500 0.6155
% 0.2400 0.3156 0.5741
% a + rand(m,n)*(b-a) 可以輸出在[a,b]之間均勻分布的隨機數組成的m行n列的矩陣。
-2 + rand(3,2) * (2 - (-2)) % 輸出在[-2,2]之間均勻分布的隨機數組成的3行2列的矩陣。
% -1.2743 0.6013
% -1.3084 0.0766
% 1.5075 0.7563
% a + rand(m,n)*(b-a)等價於unifrnd(a,b,m,n)
unifrnd(-2,2,3,2)
%% （2）代碼部分
l = 0.520; % 針的長度（任意給的）
a = 1.314; % 平行線的寬度(大於針的長度l即可)
n = 1000000; % 做n次投針試驗，n越大求出來的pi越准確
m = 0; % 記錄針與平行線相交的次數
x = rand(1, n) * a / 2 ; % 在[0, a/2]內服從均勻分布隨機產生n個數， x中每一個元素表示針的中點和最近的一條平行線的距離
phi = rand(1, n) * pi; % 在[0, pi]內服從均勻分布隨機產生n個數，phi中的每一個元素表示針和最近的一條平行線的夾角
% axis([0,pi, 0,a/2]); box on; % 畫一個坐標軸的框架，x軸位於0-pi，y軸位於0-a/2， 並打開圖形的邊框
for i=1:n % 開始循環，依次看每根針是否和直線相交
if x(i) <= l / 2 * sin(phi (i)) % 如果針和平行線相交
m = m + 1; % 那么m就要加1
% plot(phi(i), x(i), 'r.') % 模仿書上的那個圖，橫坐標為phi，縱坐標為x , 用紅色的小點進行標記
% hold on % 在原來的圖形上繼續繪制
end
end
p = m / n; % 針和平行線相交出現的頻率
mypi = (2 * l) / (a * p); % 我們根據公式計算得到的pi
disp(['蒙特卡羅方法得到pi為：', num2str(mypi)])
%% （3） 由於一次模擬的結果具有偶然性，因此我們可以重復100次后再來求一個平均的pi
result = zeros(100,1); % 初始化保存100次結果的矩陣
l = 0.520; a = 1.314;
n = 1000000;
for num = 1:100
m = 0;
x = rand(1, n) * a / 2 ;
phi = rand(1, n) * pi;
for i=1:n
if x(i) <= l / 2 * sin(phi (i))
m = m + 1;
end
end
p = m / n;
mypi = (2 * l) / (a * p);
result(num) = mypi; % 把求出來的myphi保存到結果矩陣中
end
mymeanpi = mean(result); % 計算result矩陣中保存的100次結果的均值
disp(['蒙特卡羅方法得到pi為：', num2str(mymeanpi)])
蒙特卡羅方法概述
定義
蒙特卡羅⽅法⼜稱統計模擬法，是⼀種隨機模擬⽅法，以概率和統計理論⽅法為基礎的⼀種計算⽅法，是使⽤隨機數 （或更常⻅的偽隨機數）來解決很多計算問題的⽅法。將所求解的問題同⼀定的概率模型相聯系，⽤電⼦計算機實現統計模擬或抽樣，以獲得問題的近似解。為象征性地表明這⼀⽅法的概率統計特征，故借⽤賭城蒙特卡羅命名。
提出
蒙特卡羅⽅法於20世紀40年代美國在第⼆次世界⼤戰中研制原⼦彈的“曼哈頓計划”計划的成員S.M.烏拉姆和J. 馮·諾伊曼⾸先提出。數學家馮·諾伊曼⽤馳名世界的賭城—摩納哥的Monte Carlo—來命名這種⽅法，為它蒙上了⼀ 層神秘⾊彩。在這之前，蒙特卡羅⽅法就已經存在。1777年，法國Buffon提出⽤投針實驗的⽅法求圓周率，這被認 為是蒙特卡羅⽅法的起源。
原理
由⼤數定理可知，當樣本容量⾜夠⼤時，事件的發⽣頻率即為其概率。
討論
-
蒙特卡羅是⼀種算法嗎 ?
算法（Algorithm）是指解題⽅案的准確⽽完整的描述，是⼀系列解決問題的清晰指令。蒙特卡羅准確的來說只是⼀ 種思想，或者是是⼀種⽅法。如果我們所求解的問題與概率模型有⼀定的關聯，那么我們就可以使⽤計算機多次模擬事 件發⽣，以獲得問題的近似解。從數學建模⻆度來看，⼤家千萬別認為蒙特卡羅有⼀個通⽤的代碼。每個問題對應的代 碼都是不同的，我們分析清楚題⽬后，就要⾃⼰進⾏編寫適⽤於這個題⽬的代碼。
-
蒙特卡羅與計算機仿真的關系
計算機仿真（模擬）早期稱為蒙特卡羅⽅法，是⼀⻔利⽤隨機數實驗求解隨機問題的⽅法，其主要應⽤在復雜問題的數 值模擬上。但隨着計算機的性能的提⾼以及各種新興產業的發展，⽬前計算機仿真涉及的內容要⼴得多，例如過程控制、動 畫仿真、圖像靜態模擬等都屬於計算機仿真的應⽤（例如⽤計算機模擬原⼦彈爆炸的過程、使⽤計算機⽣成特效⼤⽚等）。 在數學建模中，我們不⽤刻意的去區分兩者的區別，⼤家只需要知道如何使⽤計算機對問題進⾏模擬即可。
-
蒙特卡羅與枚舉法
枚舉法是我們中學就接觸的算法，就是把所有可能發⽣情況都考慮進去，最終計算出來⼀個確定結果。這就與蒙特卡羅 ⽅法的想法很類似，蒙特卡羅法模擬的次數越多，計算的就越准確。由於⽣活中有許多事件發⽣的結果都有⽆限種可能（例如⼀個連續分布的取值）,因此我們不可能枚舉出所有的結果，這時候就只能通過蒙特卡羅模擬，將⼀個不確定性的問題轉 化成很多個確定性問題，並得到⼀個近似解，因此蒙特卡羅算法也可以看成是枚舉法的⼀種變異。
已經學過的例⼦
第⼀期 視頻 (AOL) 第七講 : 多元回歸分析
探究內⽣性對回歸結果的影響
第⼀期 視頻 (AOL) 番外篇 : 基於熵權法對印的模型的修正
信息熵和標准差的關系
蒙特卡羅⽅法 的應⽤實例
三⻔問題
你參加⼀檔電視節⽬，節⽬組提供了ABC三扇⻔， 主持⼈告訴你，其中⼀扇⻔后邊有輛汽⻋，其它兩扇⻔ 后是空的。 假如你選擇了B⻔，這時，主持⼈打開了C⻔，讓你 看到C⻔后什么都沒有，然后問你要不要改選A⻔？
有約束的⾮線性規划問題
代碼
%% 蒙特卡羅求解有約束的非線性規划問題
% max f(x) = x1*x2*x3
% s.t.
% (1) -x1+2*x2+2*x3>=0
% (2) x1+2*x2+2*x3<=72
% (3) x2<=20 & x2>=10
% (4) x1-x2 == 10
%% （1）預備知識
% (1) format long g 可以將Matlab的計算結果顯示為一般的長數字格式（默認會保留四位小數，或使用科學計數法）
5/7
5895*514100
format long g
5/7
5895*514100
% (2)unifrnd(a,b,m,n)可以輸出在[a,b]之間均勻分布的隨機數組成的m行n列的矩陣。(等價於 a + rand(m,n)*(b-a))
unifrnd(0,5,4,3)
% 4.07361843196589 3.16179623112705 4.78753417717149
% 4.5289596853781 0.487702024997048 4.82444267599638
% 0.63493408146753 1.39249109433524 0.788065408387741
% 4.5668792806951 2.73440759602492 4.85296390880308
%% （2）代碼部分
clc,clear;
tic %計算tic和toc中間部分的代碼的運行時間
n=10000000; %生成的隨機數組數
x1=unifrnd(20,30,n,1); % 生成在[20,30]之間均勻分布的隨機數組成的n行1列的向量構成x1
x2=x1 - 10;
x3=unifrnd(-10,16,n,1); % 生成在[-10,16]之間均勻分布的隨機數組成的n行1列的向量構成x3
fmax=-inf; % 初始化函數f的最大值為負無窮（后續只要找到一個比它大的我們就對其更新）
for i=1:n
x = [x1(i), x2(i), x3(i)]; %構造x向量, 這里千萬別寫成了：x =[x1, x2, x3]
if (-x(1)+2*x(2)+2*x(3)>=0) & (x(1)+2*x(2)+2*x(3)<=72) % 判斷是否滿足條件
result = x(1)*x(2)*x(3); % 如果滿足條件就計算函數值
if result > fmax % 如果這個函數值大於我們之前計算出來的最大值
fmax = result; % 那么就更新這個函數值為新的最大值
X = x; % 並且將此時的x1 x2 x3保存到一個變量中
end
end
end
disp(strcat('蒙特卡羅模擬得到的最大值為',num2str(fmax)))
disp('最大值處x1 x2 x3的取值為：')
disp(X)
toc %計算tic和toc中間部分的代碼的運行時間
%% （3）縮小范圍重新模擬得到更加精確的取值
clc,clear;
tic %計算tic和toc中間部分的代碼的運行時間
n=10000000; %生成的隨機數組數
x1=unifrnd(22,23,n,1); % 生成在[22,23]之間均勻分布的隨機數組成的n行1列的向量構成x1
x2=x1 - 10;
x3=unifrnd(11,13,n,1); % 生成在[11,13]之間均勻分布的隨機數組成的n行1列的向量構成x3
fmax=-inf; % 初始化函數f的最大值為負無窮（后續只要找到一個比它大的我們就對其更新）
for i=1:n
x = [x1(i), x2(i), x3(i)]; %構造x向量, 這里千萬別寫成了：x =[x1, x2, x3]
if (-x(1)+2*x(2)+2*x(3)>=0) & (x(1)+2*x(2)+2*x(3)<=72) % 判斷是否滿足條件
result = x(1)*x(2)*x(3); % 如果滿足條件就計算函數值
if result > fmax % 如果這個函數值大於我們之前計算出來的最大值
fmax = result; % 那么就更新這個函數值為新的最大值
X = x; % 並且將此時的x1 x2 x3保存到一個變量中
end
end
end
disp(strcat('蒙特卡羅模擬得到的最大值為',num2str(fmax)))
disp('最大值處x1 x2 x3的取值為：')
disp(X)
toc %計算tic和toc中間部分的代碼的運行時間
書店 買書 問題 ( 0 - 1 規划)(無腦)
%% 書店買書問題的蒙特卡羅的模擬
%% （1）預備知識
% (1)unique函數: 剔除一個矩陣或者向量的重復值，並將結果按照從小到大的順序排列
% adj. 唯一的; 獨一無二的 [ju'ni:k]
unique([1 2 5; 6 8 9;2 4 6])
unique([5 6 8 8 4 1 6 2 2 4 8 4 5 6])
% (2)randi([a,b],m,n)函數可在指定區間[a,b]內隨機取出大小為m*n的整數矩陣
randi([-5,5],2,6)
%% （2）代碼求解
min_money = +Inf; % 初始化最小的花費為無窮大，后續只要找到比它小的就更新
min_result = randi([1, 6],1,5); % 初始化五本書都在哪一家書店購買，后續我們不斷對其更新
%若min_result = [5 3 6 2 3]，則解釋為：第1本書在第5家店買，第2本書在第3家店買，第3本書在第6家店買，第4本書在第2家店買，第5本書在第3家店買
n = 100000; % 蒙特卡羅模擬的次數
M = [18 39 29 48 59
24 45 23 54 44
22 45 23 53 53
28 47 17 57 47
24 42 24 47 59
27 48 20 55 53]; % m_ij 第j本書在第i家店的售價
freight = [10 15 15 10 10 15]; % 第i家店的運費
for k = 1:n % 開始循環
result = k([1, 6],1,5); % 在1-6這些整數中隨機抽取一個1*5的向量，表示這五本書分別在哪家書店購買
index = unique(result); % 在哪些商店購買了商品，因為我們等下要計算運費
money = sum(freight(index)); % 計算買書花費的運費
% 計算總花費：剛剛計算出來的運費 + 五本書的售價
for i = 1:5
money = money + M(result(i),i);
end
if money < min_money % 判斷剛剛隨機生成的這組數據的花費是否小於最小花費，如果小於的話
min_money = money % 我們更新最小的花費
min_result = result % 用這組數據更新最小花費的結果
end
end
min_money % 18+39+48+17+47+20
min_result
導彈追蹤問題
%% 蒙特卡羅用於模擬導彈追擊問題
clear;clc
%% （1）預備知識
% mod(m,n)表示求m/n的余數
mod(8,3)
mod(1000,50)
% 設置橫縱坐標的范圍並標上字符
x = 1:0.01:3;
y = x .^ 2;
plot(x,y) % 畫出x和y的圖形
axis([0 3 0 10]) % 設置橫坐標范圍為[0, 3] 縱坐標范圍為[0, 10]
pause(3) % 暫停3秒后再繼續接下來的命令
text(2,4,'清風') % 在坐標為(2,4)的點上標上字符串：清風
close % 關閉圖形窗口
%% (2) 代碼求解
% 1. 不畫追擊的示意圖
clear;clc
v=200; % 任意給定B船的速度（后期我們可以再改的）
dt=0.0000001; % 定義時間間隔
x=[0,20]; % 定義導彈和B船的橫坐標分別為x(1)和x(2)
y=[0,0]; % 定義導彈和B船的縱坐標分別為y(1)和y(2)
t=0; % 初始化導彈擊落B船的時間
d=0; % 初始化導彈飛行的距離
m=sqrt(2)/2; % 將sqrt(2)/2定義為一個常量，使后面看起來很簡潔
dd=sqrt((x(2)-x(1))^2+(y(2)-y(1))^2); % 導彈與B船的距離
while(dd>=0.001) % 只要兩者的距離足夠大，就一直循環下去。（兩者距離足夠小時表示導彈擊中，這里的臨界值要結合dt來取，否則可能導致錯過交界處的情況）
t=t+dt; % 更新導彈擊落B船的時間
d=d+3*v*dt; % 更新導彈飛行的距離
x(2)=20+t*v*m; y(2)=t*v*m; % 計算新的B船的位置 （注：m=sqrt(2)/2）
dd=sqrt((x(2)-x(1))^2+(y(2)-y(1))^2); % 更新導彈與B船的距離
tan_alpha=(y(2)-y(1))/(x(2)-x(1)); % 計算斜率，即tan(α)
cos_alpha=sqrt(1/(1+tan_alpha^2)); % sec(α)^2 = (1+tan(α)^2)
sin_alpha=sqrt(1-cos_alpha^2); % sin(α)^2 +cos(α)^2 = 1
x(1)=x(1)+3*v*dt*cos_alpha; y(1)=y(1)+3*v*dt*sin_alpha; % 計算新的導彈的位置
if d>50 % 導彈的有效射程為50個單位
disp('導彈沒有擊中B船');
break; % 退出循環
end
if d<=50 & dd<0.001 % 導彈飛行的距離小於50個單位且導彈和B船的距離小於0.001（表示擊中）
disp(['導彈飛行',num2str(d),'單位后擊中B船'])
disp(['導彈飛行的時間為',num2str(t*60),'分鍾'])
end
end
% 2. 畫追擊的示意圖
clear;clc
v=200; % 任意給定B船的速度（后期我們可以再改的）
dt=0.0000001; % 定義時間間隔
x=[0,20]; % 定義導彈和B船的橫坐標分別為x(1)和x(2)
y=[0,0]; % 定義導彈和B船的縱坐標分別為y(1)和y(2)
t=0; % 初始化導彈擊落B船的時間
d=0; % 初始化導彈飛行的距離
m=sqrt(2)/2; % 將sqrt(2)/2定義為一個常量，使后面看起來很簡潔
dd=sqrt((x(2)-x(1))^2+(y(2)-y(1))^2); % 導彈與B船的距離
for i=1:2
plot(x(i),y(i),'.k','MarkerSize',1); % 畫出導彈和B船所在的坐標，點的大小為1，顏色為黑色(k)，用小點表示
grid on; % 打開網格線
hold on; % 不關閉圖形，繼續畫圖
end
axis([0 30 0 10]) % 固定x軸的范圍為0-30 固定y軸的范圍為0-10
k = 0; % 引入一個變量 為了控制畫圖的速度（因為Matlab中畫圖的速度超級慢）
while(dd>=0.001) % 只要兩者的距離足夠大，就一直循環下去。（兩者距離足夠小時表示導彈擊中，這里的臨界值要結合dt來取，否則可能導致錯過交界處的情況）
t=t+dt; % 更新導彈擊落B船的時間
d=d+3*v*dt; % 更新導彈飛行的距離
x(2)=20+t*v*m; y(2)=t*v*m; % 計算新的B船的位置 （注：m=sqrt(2)/2）
dd=sqrt((x(2)-x(1))^2+(y(2)-y(1))^2); % 更新導彈與B船的距離
tan_alpha=(y(2)-y(1))/(x(2)-x(1)); % 計算斜率，即tan(α)
cos_alpha=sqrt(1/(1+tan_alpha^2)); % 利用公式：sec(α)^2 = (1+tan(α)^2) 計算出cos(α)
sin_alpha=sqrt(1-cos_alpha^2); % 利用公式： sin(α)^2 +cos(α)^2 = 1 計算出sin(α)
x(1)=x(1)+3*v*dt*cos_alpha; y(1)=y(1)+3*v*dt*sin_alpha; % 計算新的導彈的位置
k = k +1 ;
if mod(k,500) == 0 % 每刷新500次時間就畫出下一個導彈和B船所在的坐標 mod(m,n)表示求m/n的余數
for i=1:2
plot(x(i),y(i),'.k','MarkerSize',1);
hold on; % 不關閉圖形，繼續畫圖
end
pause(0.001); % 暫停0.001s后再繼續下面的操作
end
if d>50 % 導彈的有效射程為50個單位
disp('導彈沒有擊中B船');
break; % 退出循環
end
if d<=50 & dd<0.001 % 導彈飛行的距離小於50個單位且導彈和B船的距離小於0.001（表示擊中）
disp(['導彈飛行',num2str(d),'個單位后擊中B船'])
disp(['導彈飛行的時間為',num2str(t*60),'分鍾'])
end
end
作業
作業1
%% 蒙特卡羅求解非線性規划問題
% min f(x) =2*(x1^2)+x2^2-x1*x2-8*x1-3*x2
% s.t.
% (1) 3*x1+x2>9
% (2) x1+2*x2<16
% (3) x1>0 & x2>0
%% （1）初次尋找最小值的代碼
clc,clear;
format long g %可以將Matlab的計算結果顯示為一般的長數字格式（默認會保留四位小數，或使用科學計數法）
tic %計算tic和toc中間部分的代碼的運行時間
n=10000000; %生成的隨機數組數
x1=unifrnd(0,16,n,1); % 生成在[0,16]之間均勻分布的隨機數組成的n行1列的向量構成x1
x2=unifrnd(0,8,n,1); % 生成在[0,8]之間均勻分布的隨機數組成的n行1列的向量構成x2
fmin=+inf; % 初始化函數f的最小值為正無窮（后續只要找到一個比它小的我們就對其更新）
for i=1:n
x = [x1(i), x2(i)]; %構造x向量, 這里千萬別寫成了：x =[x1, x2]
if (3*x(1)+x(2)>9) & (x(1)+2*x(2)<16) % 判斷是否滿足條件
result = 2*(x(1)^2)+x(2)^2-x(1)*x(2)-8*x(1)-3*x(2); % 如果滿足條件就計算函數值
if result < fmin % 如果這個函數值小於我們之前計算出來的最小值
fmin = result; % 那么就更新這個函數值為新的最小值
X = x; % 並且將此時的x1 x2 保存到相應的變量中
end
end
end
disp(strcat('蒙特卡羅模擬得到的最小值為',num2str(fmin)))
disp('最小值處x1 x2的取值為：')
disp(X)
toc %計算tic和toc中間部分的代碼的運行時間
%% （2）縮小范圍重新模擬得到更加精確的取值
clc,clear;
tic %計算tic和toc中間部分的代碼的運行時間
n=10000000; %生成的隨機數組數
x1=unifrnd(2,3,n,1); % 生成在[2,3]之間均勻分布的隨機數組成的n行1列的向量構成x1
x2=unifrnd(2,3,n,1); % 生成在[2,3]之間均勻分布的隨機數組成的n行1列的向量構成x2
fmin=+inf; % 初始化函數f的最小值為正無窮（后續只要找到一個比它小的我們就對其更新）
for i=1:n
x = [x1(i), x2(i)]; %構造x向量, 這里千萬別寫成了：x =[x1, x2]
if (3*x(1)+x(2)>9) & (x(1)+2*x(2)<16) % 判斷是否滿足條件
result = 2*(x(1)^2)+x(2)^2-x(1)*x(2)-8*x(1)-3*x(2); % 如果滿足條件就計算函數值
if result < fmin % 如果這個函數值小於我們之前計算出來的最小值
fmin = result; % 那么就更新這個函數值為新的最小值
X = x; % 並且將此時的x1 x2 保存到相應的變量中
end
end
end
disp(strcat('蒙特卡羅模擬得到的最小值為',num2str(fmin)))
disp('最小值處x1 x2的取值為：')
disp(X)
toc %計算tic和toc中間部分的代碼的運行時間
作業2
某設備上安裝有四只型號規格完全相同的電子管，已知電子管壽命（假定為整數）為1000h~2000h之間的均勻分布.當電子管損壞時有兩種維修方案，一是每次更換損壞的那一只；二是當其中一只損壞時四只同時更換.已知更換時間為換一只時需1h，4只同時換為2h.更換時機器因停止運轉每小時的損失為20元，又每只電子管價格10元，試用模擬方法決定哪一個方案經濟合理？
%% 選擇決策方案的模擬
% 某設備上安裝有四只型號規格完全相同的電子管，已知電子管壽命為1000--2000小時之間的均勻分布(假定為整數)。
% 當電子管損壞時有兩種維修方案，一是每次更換損壞的那一只；二是當其中一只損壞時四只同時更換。
% 已知更換時間為換一只時需1小時，4只同時換為2小時。
% 更換時機器因停止運轉每小時的損失為20元，又每只電子管價格10元，
% 試用模擬方法決定哪一個方案經濟合理？
%% （1）預備知識
% randi([a,b],m,n) 隨機生成m*n的矩陣，矩陣中的每個元素都是[a,b]中的隨機整數
randi([1, 5],3,2)
randi([1, 5]) % 不寫m*n代表只生成1個隨機數
% find函數的用法
% find函數的用法在第一期視頻：層次分析法那一節講過，我們當時找最大特征值的位置
a = [2 3 5 1 7 5];
find(a) % 找到a中所有非0元素的位置
find(a == 5) % 找到a中等於5的元素的位置
find(a == 5,1) % 找到a中第一個等於5的元素的位置
find(a == min(a)) % 找到a中最小元素的位置
%% （2）代碼部分
clear;clc
T = 100000000; % T表示模擬的總時間（單位為小時）
t = 0; % 初始化當前時刻為0小時
c1 = 0; c2 = 0; % 初始化兩種方案的總花費都為0
%% 方案一
life = randi([1000,2000],1,4); % 隨機生成四個電子管的壽命，假設為整數
while t < T % 只要現在的時刻沒有超過總時刻，就不斷循環下去
result = min(life); % 找出壽命最短的那一個電子管的壽命
t = t+result+1; % 現在的時間更改到有電子管損壞的時刻（加上1表示更換電子管需要花費的時間）
c1 = c1 + 20 * 1 +10; % 更新方案一的花費
k = find(life == result,1); % 找到哪一個電子管是壞的
life = life - result -1; % 更新所有電子管的壽命
life(k) = randi([1000,2000]); % 把壞掉的那個電子管的壽命重置
end
%% 方案二
t = 0; % 初始化當前時刻為0小時
while t < T % 只要現在的時刻沒有超過總時刻，就不斷循環下去
life = randi([1000,2000],1,4); % 隨機生成四個電子管的壽命，假設為整數
result = min(life); % 找出壽命最小的那一個電子管的壽命
t = t+result+2; % 現在的時間更改到有電子管損壞的時刻（加上2表示更換所有電子管需要花費的時間）
c2 =c2 + 20 * 2 +40; % 更新方案二的花費
end
%% 兩種方案的花費
c1
c2