汽油与消费需求问题的MonteCarlo求解方法

一.问题的提出
二.问题的分析
三.代码的实现

一.问题的提出

你受雇于一家加油站连锁店当顾问，你现在的任务是要确定每隔多长时间把多少汽油运送到各个加油站。每次运送汽油都要支付费用ddd，其是于运送无关的附加费用。加油站都在告诉公路旁，所以需求量可以合理的假设为常数，决定费用的其他因素有:

存储中冻结的资金
分期偿还的设备费用
保险费
税费
安全检测费

假定在短期运营中每个加油站汽油的需求和价格都是常数。只要加油站没有短缺，就可以得到不变的总收入。因此如果我们想要最大化总利润，那么就必须最小化总费用。提出以下的问题:
在每个加油站存储足够汽油的同时，最小化每天的平均费用。

二.问题的分析

首先，我们根据常识可以建立以下的模型：
日平均费用=f(存储费用，运费，需求率)日平均费用 = f(存储费用，运费，需求率) 日平均费用=f(存储费用，运费，需求率)
在这里我们可以合理的假设，在问题设计的范围内，单位存储量费用是常数。且日需求量q‾\overline qq是定值。

理论上来说，我们分析该日平均费用模型，可以得到一个最优的运送时间间隔和最优运量：
T∗=2dsrQ=rT∗T^* = \sqrt{\frac{2d}{sr}}\\ Q = rT^* T∗=sr2dQ=rT∗
其中，T∗:T^*:T∗:最佳时间运送时间间隔(天)，Q∗:Q^*:Q∗:汽油最佳运量(加仑)，d:d:d:每次运送的费用(元)，s:s:s:每加仑汽油每天的存储费(元/天)，r:r:r:汽油每日的需求量(加仑/天)。

上面的式子就是理论考虑的过程，我们暂不考虑，只是单纯的对这个模型进行以下仿真。

正如我们前面所讨论的，montecarlo方法最重要的一步就是建立实际模型中的随机因素和matlab生成随机数之间的一个映射关系。很明显，一个加油站每天特定的需求量rrr很明显应该是一个随机因素，如果我们通过过去1000天的研究调查发现其每天的需求量和出现的天数的表格如下:

日需求量rrr(加仑/天)	出现天数	出现概率
1000~1099	10	0.01
1100~1199	20	0.02
1200~1299	50	0.05
1300~1399	120	0.12
1400~1499	200	0.20
1500~1599	270	0.27
1600~1699	180	0.18
1700~1799	80	0.08
1800~1899	40	0.04
1900~1999	30	0.03

如果我们有一个在[0,1][0,1][0,1]上的随机数我们可以控制其生成的随机数若落在在不同长度的子区间上，来控制其不同的概率值。比如我们要模拟产生1100~1199之间的一个随机的日需求量，那么此时我们可以产生一个[0,1][0,1][0,1]之内的随机数，如果其正好落在[0.01,0.02][0.01 ,0.02][0.01,0.02]区间内，则说明的确模拟产生了一个1100~1199之内的一个日需求量，但是到底是多少，我们就需要在这个区间内进行线性样条插值(或者三次样条插值)。目的就是为了使得我们能用[0,1][0,1][0,1]区间之内的一个随机数去成功映射到一个日需求量上:
f−1:rand(0,1)→rf^{-1}:rand(0,1)\rightarrow r f−1:rand(0,1)→r
那么我们这个映射f:r→rand(0,1)f:r\rightarrow rand(0,1)f:r→rand(0,1)用样条插值建立出来如下：

需求区间	线性样条
1000~1050	s(q)=0.0002q−0.2s(q) = 0.0002q-0.2s(q)=0.0002q−0.2
1050~1150	s(q)=0.0002q−0.2s(q) = 0.0002q-0.2s(q)=0.0002q−0.2
1150~1250	s(q)=0.0005q−0.545s(q) = 0.0005q-0.545s(q)=0.0005q−0.545
1250~1350	s(q)=0.0012q−1.42s(q) = 0.0012q-1.42s(q)=0.0012q−1.42
1350~1450	s(q)=0.002q−2.5s(q) = 0.002q-2.5s(q)=0.002q−2.5
1450~1550	s(q)=0.0027q−3.515s(q) = 0.0027q-3.515s(q)=0.0027q−3.515
1550~1650	s(q)=0.0018q−2.12s(q) = 0.0018q-2.12s(q)=0.0018q−2.12
1650~1750	s(q)=0.0008q−0.47s(q) = 0.0008q-0.47s(q)=0.0008q−0.47
1750~1850	s(q)=0.0004q+0.23s(q) = 0.0004q+0.23s(q)=0.0004q+0.23
1850~1950	s(q)=0.0002q+0.6s(q) = 0.0002q+0.6s(q)=0.0002q+0.6

这个式子求出来了之后，我们要根据随机数求rrr只需要将上面的式子求逆运算即可。

现在，让我们回归正题，考虑如何用我们上面建立的rrr经验数据进行monteCarlo模拟：

符号变量	含义
QQQ	汽油运量(加仑)
TTT	运送的时间间隔(天)
III	当前的存储量(加仑)
ddd	每次运送的费用(元)
sss	每加仑汽油每天的存储费用(元/天*加仑)
CCC	总费用(元)
ccc	每天的平均费用(元/天)
NNN	模拟运行的天数（天）
KKK	模拟中剩下的天数(天)
ζi\zeta_iζi	[0,1][0,1][0,1]区间的随机数
qiq_iqi	日需求量(加仑/天)
FlagFlagFlag	终止计算的标志

INPUTS:Q,T,d,s,NINPUTS:Q,T,d,s,NINPUTS:Q,T,d,s,N

OUTPUT:cOUTPUT:cOUTPUT:c

STEP1: 初始化K=N,I=0,C=0,Flag=0K=N,I=0,C = 0,Flag = 0K=N,I=0,C=0,Flag=0

STEP2: 开始一次运送后开始下个存储期I=I+Q,C=C+dI = I+Q,C = C+dI=I+Q,C=C+d

STEP3: 判断如果T=KT=KT=K,则Flag=1Flag = 1Flag=1

STEP4: 接下来求这个存储期内的费用：i=1,2,3...Ti=1,2,3...Ti=1,2,3...T

STEP5: 根据ζi\zeta_iζi计算qiq_iqi，并且修正当前的存储量I=I−qiI = I - q_iI=I−qi

STEP6:判断如果I≤0I\leq0I≤0,表示当前的存储量全部用完，就让I=0I = 0I=0,转到STEP8,否则转到STEP7

STEP7: 计算日存储和总费用：C=C+IsC = C + IsC=C+Is

STEP8: 天数K=K−1K = K-1K=K−1

STEP9：如果Flag=0Flag = 0Flag=0转到STEP2,否则，转下一步。

STEP10: 计算c=CNc = \frac{C}{N}c=NC并且输出。

三.代码的实现

我们设N=36500,I=0,C=0,Q=3000,T=3,s=2,d=1000N = 36500,I = 0,C = 0,Q = 3000,T = 3,s = 2,d = 1000N=36500,I=0,C=0,Q=3000,T=3,s=2,d=1000。

function zhongzi
N = 36500;
I = 0;C = 0;
Q = 3000;
c = zeros(1,N);
T = 3;d = 1000;s = 2;
for i = 1:T:NI = I + Q;C = C + d;for j = 1:Tq = randTrans(unifrnd(0,1));I = I - q;if I<=0I = 0;break;endC = C + I*s;c(i+j) = C/(i+j);end
end
c(1) = c(2);
for j = 2:Nif c(j) == 0c(j) = c(j-1);end
end
[m,N] = size(c);
hold on;
plot(1:N,c,'b');
plot(1:N,c(N)*ones(1,N),'r');
xlabel('天数');
ylabel('日平均费用');
title('仿真收敛曲线');
legend('实际值','最终值');
disp(['日平均费用大概收敛于：',num2str(c(N))]);
endfunction q = randTrans(x)
k = [5000 5000 2000 833.33 500 370.37 555.55 1250 2500 5000];
x0 = [0.2 0.2 0.545 1.42 2.5 3.515 2.12 0.47 -0.23 -0.6];
if 0<=x<0.01q = (x+x0(1))*k(1);
elseif 0.01<=x<0.03q = (x+x0(2))*k(2);
elseif 0.03<=x<0.08    q = (x+x0(3))*k(3);
elseif 0.08<=x<0.20q = (x+x0(4))*k(4);
elseif 0.20<=x<0.40q = (x+x0(5))*k(5);
elseif 0.40<=x<0.67q = (x+x0(6))*k(6)
elseif 0.67<=x<0.85q = (x+x0(7))*k(7);
elseif 0.85<=x<0.93q = (x+x0(8))*k(8);
elseif 0.93<=x<0.97q = (x+x0(9))*k(9);
elseq = (x+x0(10))*k(10);
end
end

运行之后得到的结果是:

日平均费用大概收敛于：1435.8289