###### tags: `Python筆記`
# 【Python】決策科學作業4
## 第一題
```python
%pip install scipy --upgrade
import math
import matplotlib.pyplot as plt
import scipy.stats as stats
import scipy.optimize as opt
import numpy as np
import pandas as pd
%pip install statsmodels --upgrade
import statsmodels.api as sm
%pip install pymoo
import pymoo as pymoo
```
```python
#homework demand function:
arrD = [0,1,2,3,4,5,6,7,8,9,10]
def Demand(T):
d = np.random.choice(arrD , T , p=[0.01,0.02,0.04,0.06,0.09,0.14,0.18,0.22,0.16,0.06,0.02] )
return d
#Holding cost per item per day
h_val=0.6
unit_price=400
unit_cost=250
#Profit margin per TV sold
unit_profit=unit_price-unit_cost
#Cost per delivery from supplier to the dealer
shipping_cost=500
def InvModelsS(x):
#Set the seed for reproducibility
np.random.seed(5566)
#you have two decision variables
s=x[0]
S=x[1]
#The number of simulation runs
M=5000
#The number of days
T=360
#The experiment will be repeated for
#m in 1:M repetitions & t in 1:T duration
#Vector initialization
#SL: service level
SL=np.empty(M)
daily_profit=np.empty(M)
for m in range(M):
inv_onhand=np.empty(T+1)
inv_onhand[0]=S
sales=np.zeros(T)
orderQ=np.zeros(T)
#Keep track of inventory on-order but not arrived yet
inv_onorder=np.array([])
inv_onorder_stamp=np.array([])
#The number of deliveries from the supplier
delivery=0
#The total amount of lost sales
loss=0
#The number of stockouts in a repetition
stockout=0
#Simulate Poisson demand for T periods
# d=np.random.poisson(lambda_est,T)
d = Demand(T)
arrival=0
for t in range(T):
if t==0:
#Selling process
sales[t]=min(inv_onhand[t],d[t])
#Check stockout
if (d[t]>sales[t]):
stockout=stockout+1
loss=loss+(d[t]-sales[t])
#Compute inventory position
inv_position=inv_onhand[0]
#Ordering mechanism of (s, S)
if inv_position<=s:
#Compute order quantity
#This hw need to order the fixed quantity of Q
Q=S
#Update inventory on-order
inv_onorder=np.append(inv_onorder,Q)
#Simulate stochastic lead time
#change the lead time prob in homework
time_to_arrive=t+np.random.choice([3,4,5],size=1,p=[0.2,0.6,0.2])+1
#Update time stamp for invenory on-order
inv_onorder_stamp=np.append(inv_onorder_stamp,time_to_arrive)
#Record order quantity
orderQ[t]=Q
#Update inventory on-hand in the end of each day
inv_onhand[t+1]=inv_onhand[t]-sales[t]
###End t==0
if t>0:
#Check if any inventory on-order should arrive
if np.any(inv_onorder_stamp==t):
#Update the number of deliveries
delivery=delivery+1
#Compute the total of arrived inventories
index=np.where(inv_onorder_stamp==t)
arrival=np.sum(inv_onorder[index])
#Update inventory on-hand before starting the day
inv_onhand[t]=inv_onhand[t]+arrival
#Removed those just arrived from inventory on-order
inv_onorder=np.delete(inv_onorder,index)
inv_onorder_stamp=np.delete(inv_onorder_stamp,index)
#Record sales
sales[t]=min(inv_onhand[t],d[t])
#Check if any stockout takes place
if d[t]>sales[t]:
stockout=stockout+1
loss=loss+d[t]-sales[t]
#End if
#Update inventory position
inv_position=inv_onhand[t]+np.sum(inv_onorder)
#Ordering mechanism of (s, S) inventory policy
if (inv_position<=s):
#This hw need to order the fixed quantity of Q
Q=S
inv_onorder=np.append(inv_onorder,Q)
time_to_arrive=t+np.random.choice([3,4,5],size=1,p=[0.2,0.6,0.2])+1
inv_onorder_stamp=np.append(inv_onorder_stamp,time_to_arrive)
orderQ[t]=Q
#Update inventory on-hand in the end of the day
inv_onhand[t+1]=inv_onhand[t]-sales[t]
###End t>0
#print("t=",t,"d=",d[t],"arrival=",arrival,"; inv_onhand=",inv_onhand[t],"; order(t)=",orderQ[t])
###End t loop
#Calculate the service level (i.e., in-stock probability)
SL[m]=(T-stockout)/t
#Calculate the average daily profit
daily_profit[m]=(unit_profit*np.sum(sales)-np.sum(h_val*inv_onhand)-shipping_cost*delivery-unit_profit*loss)/T
#Print simulation progrss
#if (np.remainder(m,100)==0):
# print("m=",m)
return -1*np.mean(daily_profit)
```
```python
##Use Pymoo for numerical optimization
##Define the optimization problem first
from pymoo.problems.functional import FunctionalProblem
problemInv=FunctionalProblem(n_var=2,objs=InvModelsS,
xl=np.array([20,100]),xu=np.array([90,200]))
```
```python
##Hooke-Jeeves Search
from pymoo.algorithms.soo.nonconvex.pattern import PatternSearch
from pymoo.optimize import minimize
resHK=minimize(problemInv,algorithm=PatternSearch(),verbose=True)
print("Best solution found: \nX = %s\nF = %s" % (resHK.X, resHK.F))
```
```
=================================================
n_gen | n_eval | f_avg | f_min
=================================================
1 | 20 | -7.987321E+02 | -8.268439E+02
2 | 24 | -8.269384E+02 | -8.270328E+02
3 | 29 | -8.270328E+02 | -8.270328E+02
4 | 33 | -8.270328E+02 | -8.270328E+02
5 | 37 | -8.270452E+02 | -8.270576E+02
6 | 42 | -8.270576E+02 | -8.270576E+02
7 | 46 | -8.270638E+02 | -8.270700E+02
8 | 51 | -8.271015E+02 | -8.271330E+02
9 | 56 | -8.271330E+02 | -8.271330E+02
10 | 60 | -8.271330E+02 | -8.271330E+02
11 | 64 | -8.271330E+02 | -8.271330E+02
12 | 68 | -8.271330E+02 | -8.271330E+02
13 | 72 | -8.271330E+02 | -8.271330E+02
14 | 76 | -8.271330E+02 | -8.271330E+02
15 | 80 | -8.271330E+02 | -8.271330E+02
16 | 84 | -8.271330E+02 | -8.271330E+02
17 | 88 | -8.271330E+02 | -8.271330E+02
18 | 92 | -8.271330E+02 | -8.271330E+02
```
**Best solution found:**
X = [ 45.23958658 100. ]
F = [-827.13301022]
```python
##Nelder-Mead Search
from pymoo.algorithms.soo.nonconvex.nelder import NelderMead
resNM=minimize(problemInv,algorithm=NelderMead(),verbose=True)
print("Best solution found: \nX = %s\nF = %s" % (resNM.X, resNM.F))
```
```
=================================================
n_gen | n_eval | f_avg | f_min
=================================================
1 | 20 | -7.971971E+02 | -8.210567E+02
2 | 22 | -8.202593E+02 | -8.210567E+02
3 | 23 | -8.198597E+02 | -8.210567E+02
4 | 25 | -8.210112E+02 | -8.224950E+02
5 | 27 | -8.232974E+02 | -8.263406E+02
6 | 29 | -8.252921E+02 | -8.270408E+02
7 | 31 | -8.261500E+02 | -8.270408E+02
8 | 33 | -8.267295E+02 | -8.270408E+02
9 | 35 | -8.268829E+02 | -8.270408E+02
10 | 37 | -8.269563E+02 | -8.270408E+02
11 | 41 | -8.269082E+02 | -8.270408E+02
12 | 43 | -8.269913E+02 | -8.270568E+02
13 | 45 | -8.270613E+02 | -8.270862E+02
14 | 47 | -8.270701E+02 | -8.270862E+02
15 | 51 | -8.270680E+02 | -8.270944E+02
16 | 53 | -8.270919E+02 | -8.270952E+02
17 | 55 | -8.271026E+02 | -8.271181E+02
18 | 57 | -8.271237E+02 | -8.271579E+02
19 | 61 | -8.271493E+02 | -8.271579E+02
20 | 65 | -8.270960E+02 | -8.271740E+02
21 | 67 | -8.271198E+02 | -8.271740E+02
22 | 69 | -8.271329E+02 | -8.271740E+02
23 | 71 | -8.271613E+02 | -8.271740E+02
24 | 73 | -8.271632E+02 | -8.271740E+02
25 | 75 | -8.271643E+02 | -8.271740E+02
26 | 77 | -8.271770E+02 | -8.271961E+02
27 | 81 | -8.271831E+02 | -8.271961E+02
28 | 83 | -8.271894E+02 | -8.271961E+02
29 | 85 | -8.271956E+02 | -8.271961E+02
30 | 87 | -8.271970E+02 | -8.271987E+02
31 | 89 | -8.271979E+02 | -8.271988E+02
32 | 91 | -8.271983E+02 | -8.271988E+02
33 | 93 | -8.271992E+02 | -8.272000E+02
34 | 95 | -8.271993E+02 | -8.272000E+02
35 | 97 | -8.272000E+02 | -8.272010E+02
36 | 99 | -8.272010E+02 | -8.272020E+02
37 | 101 | -8.272025E+02 | -8.272045E+02
38 | 103 | -8.272039E+02 | -8.272054E+02
39 | 105 | -8.272059E+02 | -8.272079E+02
40 | 107 | -8.272074E+02 | -8.272088E+02
41 | 109 | -8.272079E+02 | -8.272088E+02
42 | 111 | -8.272081E+02 | -8.272088E+02
43 | 113 | -8.272086E+02 | -8.272091E+02
44 | 115 | -8.272093E+02 | -8.272100E+02
45 | 117 | -8.272095E+02 | -8.272100E+02
46 | 119 | -8.272095E+02 | -8.272100E+02
47 | 121 | -8.272096E+02 | -8.272100E+02
48 | 123 | -8.272097E+02 | -8.272100E+02
49 | 125 | -8.272099E+02 | -8.272101E+02
50 | 127 | -8.272100E+02 | -8.272101E+02
51 | 129 | -8.272100E+02 | -8.272101E+02
52 | 131 | -8.272101E+02 | -8.272102E+02
53 | 133 | -8.272102E+02 | -8.272102E+02
54 | 135 | -8.272102E+02 | -8.272103E+02
55 | 137 | -8.272102E+02 | -8.272103E+02
56 | 139 | -8.272103E+02 | -8.272103E+02
57 | 141 | -8.272103E+02 | -8.272103E+02
58 | 143 | -8.272103E+02 | -8.272103E+02
59 | 145 | -8.272103E+02 | -8.272103E+02
60 | 147 | -8.272103E+02 | -8.272103E+02
61 | 149 | -8.272103E+02 | -8.272103E+02
62 | 151 | -8.272103E+02 | -8.272103E+02
63 | 153 | -8.272103E+02 | -8.272103E+02
64 | 155 | -8.272103E+02 | -8.272103E+02
65 | 157 | -8.272103E+02 | -8.272103E+02
66 | 159 | -8.272103E+02 | -8.272103E+02
67 | 161 | -8.272103E+02 | -8.272103E+02
68 | 163 | -8.272103E+02 | -8.272103E+02
69 | 165 | -8.272103E+02 | -8.272103E+02
70 | 167 | -8.272103E+02 | -8.272103E+02
71 | 169 | -8.272103E+02 | -8.272103E+02
72 | 171 | -8.272103E+02 | -8.272103E+02
73 | 173 | -8.272103E+02 | -8.272103E+02
74 | 175 | -8.272103E+02 | -8.272103E+02
75 | 177 | -8.272103E+02 | -8.272103E+02
76 | 179 | -8.272103E+02 | -8.272103E+02
77 | 181 | -8.272103E+02 | -8.272103E+02
78 | 183 | -8.272103E+02 | -8.272103E+02
79 | 185 | -8.272103E+02 | -8.272103E+02
80 | 187 | -8.272103E+02 | -8.272103E+02
81 | 189 | -8.272103E+02 | -8.272103E+02
82 | 191 | -8.272103E+02 | -8.272103E+02
83 | 193 | -8.272103E+02 | -8.272103E+02
84 | 195 | -8.272103E+02 | -8.272103E+02
85 | 197 | -8.272103E+02 | -8.272103E+02
86 | 199 | -8.272103E+02 | -8.272103E+02
87 | 201 | -8.272103E+02 | -8.272103E+02
88 | 203 | -8.272103E+02 | -8.272103E+02
89 | 205 | -8.272103E+02 | -8.272103E+02
90 | 207 | -8.272103E+02 | -8.272103E+02
91 | 209 | -8.272103E+02 | -8.272103E+02
92 | 211 | -8.272103E+02 | -8.272103E+02
```
**Best solution found:**
X = [ 45.44456331 103.77779118]
F = [-827.21033316]
Sign in with Wallet
Connect another wallet