Packages Used: No new packages are used in this notebook.
c = [[10, 20, 30], 40, [50, 60]]
3-element Vector{Any}:
[10, 20, 30]
40
[50, 60]
The individual top-level elements in the array are each an array that can be accessed:
c[1]
3-element Vector{Int64}:
10
20
30
To access an element in the second-level array:
c[1][1]
10
To add an element to the end of a second-level array:
push!(c[1],35)
4-element Vector{Int64}:
10
20
30
35
To add an element to the end of the top-level array:
push!(c,[1:2:10;])
4-element Vector{Any}:
[10, 20, 30, 35]
40
[50, 60]
[1, 3, 5, 7, 9]
Arrays can be used to store text strings:
t = ["Miami","Detroit","Boston"]
3-element Vector{String}:
"Miami"
"Detroit"
"Boston"
Many functions can accept arrays of strings as input:
t = sort(t)
3-element Vector{String}:
"Boston"
"Detroit"
"Miami"
idxDetroit = findall("Detroit" .== t)
1-element Vector{Int64}:
2
If not sure of capitalization, can convert to lowercase before string matching:
isDetroit = "detroit" .== lowercase.(t)
3-element BitVector: 0 1 0
lowercase.(t)
3-element Vector{String}:
"boston"
"detroit"
"miami"
A tuple is similar to an array, but differs in that, it is defined using ( ) instead of [ ] and, once defined, it cannot be changed and is immutable. A tuple is used to display multiple outputs:
a = 1; b = 2; c = 3;
a, b, c
(1, 2, 3)
a,b,c = 1,2,3
(1, 2, 3)
a
1
A tuple can be used to store all of the input arguments for a function:
abc = (1, 2, 3)
(1, 2, 3)
abc[1]
1
Define a function of three inputs:
f(x,y,z) = x + y * z
f(a, b, c)
7
The single-input tuple abc can be used instead to pass to "splat" the inputs a, b, c to the function by using ... following the tuple:
f(abc...)
7
Based on the name Operations Research, it is not surprising that the techniques of OR are used to analyze the operations of a variety of different systems, including production systems that are involved in manufacturing goods or providing services; for example, a hospital can be thought of as a system that produces healthcare services. In most cases, the OR analysis will involve developing a model of the system and then somehow optimizing the model in the hope that these results can be transferred into improved operation of the system. Nonlinear optimization is the OR technique that has the widest application in operations analysis.
A testing center has a device that is able to use ultrasonic waves transmitted within a special solution to detect defects in parts like aircraft engines. The part is placed in a tank that, in order to operate correctly, is filled with just enough solution to submerge the part; the solution is then drained from the tank and disposed of since it can become contaminated during testing. Currently, test results are being delayed because it is time-consuming to fill the tank with the solution. In order to speed up testing, a suggestion has been made to consider pre-filling the tank and then either adding or draining the solution as needed to submerge a part. The following demand data is available regarding the cubic volume of solution that was needed to test recent parts: 237, 214, 161, 146, 331, 159, 423, 332, 139, 327, 152, 98, and 116. Determine the cubic volume of solution that should be pre-filled in the tank to minimize the time required to submerge the part, given that the tank can be filled or drained at the same rate.
It is a good idea to look at the data to see if there is a pattern to the data and check to see if there are any anomalous, outlier data that could be erroneous. If anything is found, it might be possible to ask whoever provided the data about the issue before doing any further analysis. If some data is in error, then it can be deleted; outlier data should not automatically be eliminated unless it is a real error.
It is easy to display 1-D, a little harder but doable for 2-D, and difficult for 3+ degree data, where the best you can do usually is look at different 1- and 2-D cross-sections of the day. For this data, since there is no time element, a histogram is used instead of a scatter or line plot.
using Plots
d = [237, 214, 161, 146, 331, 159, 423, 332, 139, 327, 152, 98, 116]
histogram(d, legend=false, xlabel="Demand", ylabel="Count")
Not seeing anything unusual in the histogram, the next step can be to develop a model that captures the "cost" of different pre-fill levels, where the cost, in this case, can be the time required to either fill or drain the tank from the pre-fill level, where it is assumed that, for a given volume of solution, the fill and drain times are the same since their rates are the same.
$ \quad \mbox{Time: } t(q,d_i) = \bigl| d_i - q \bigr|$
ft(q,di) = abs(di - q) # Time for single part
fTT(q,d) = sum([ft(q,di) for di in d]) # Total time
plot(q -> fTT(q,d), 0, maximum(d), legend=false, xlabel="Pre-Fill Level", ylabel="Total Time")
using Optim
qᵒ = optimize(q -> fTT(q,d), 0, maximum(d)).minimizer # (qᵒ ⟹ q\^o<TAB>)
display(scatter!([qᵒ],[fTT(qᵒ,d)]))
qᵒ, fTT(qᵒ,d)
(161.0000005146973, 1054.0000005146974)
Since the solution is only as good as the model and data used to generate it, it is helpful to try to determine how robust the solution is with respect to the data that was used and the assumptions behind the model.
Starting with data robustness, one thing that can be checked is what if there had been an extreme outlier in the data. For example, if the 423 demand was instead 1423, would the solution change significantly?
d2 = copy(d) # Note: Need to use `copy` instead of d2 = d; otherwise d2 just
d2[d2 .== 423] .= 1423 # points to d and any change to d2 will also change d
@show d d2;
d = [237, 214, 161, 146, 331, 159, 423, 332, 139, 327, 152, 98, 116] d2 = [237, 214, 161, 146, 331, 159, 1423, 332, 139, 327, 152, 98, 116]
plot!(q -> fTT(q,d2), 0, maximum(d2))
qᵒ = optimize(q -> fTT(q,d2), 0, maximum(d2)).minimizer
display(scatter!([qᵒ],[fTT(qᵒ,d2)]))
qᵒ, fTT(qᵒ,d2)
(160.999999114814, 2054.000000885186)
Although, as expected, the total time increases by 1000, it might seem surprising that the optimal pre-fill level remains unchanged. The reason for this is, for this particular objective function, the optimal solution is the value of the median demand point. Any unit increase in the solution would increase by one unit the times of all of the demands less and decrease by one unit the times of all demands greater than the solution. Since the unit cost of filling and draining are the same and, being the median point, there are equal numbers of points less than and greater than the solution, the optimal solution remains unchanged.
using Statistics
median(d), median(d2)
(161.0, 161.0)
The assumption in the model that is most questionable is assuming that the fill and train rates are the same; this would only be true if the diameters of the fill and drain tubes are the same along with the flow rates being the same. If, after further investigation, it was found that the drain rate was twice the fill rate, then the model can be modified to incorporate this. Since the "cost" of filling and raining is now unbalanced, with filling taking twice as long as draining, the median point will not, in general, correspond to the optimal solution.
$ \quad \mbox{Pre-fill cost: } c(q,d_i) = \begin{cases} 2(d_i - q), & \mbox{if } d_i > q \\ q - d_i, & \mbox{otherwise } \end{cases}$
fc(q,di) = di > q ? 2(di - q) : q - di # Time for single part
fTC(q,d) = sum([fc(q,di) for di in d]) # Total time
qᵒ = optimize(q -> fTC(q,d), 0, maximum(d)).minimizer
q2ᵒ = optimize(q -> fTC(q,d2), 0, maximum(d2)).minimizer
qᵒ,q2ᵒ
(237.00000215131104, 237.00000093451362)
As expected, since it is more costly to fill than to drain, the optimal pre-fill level has increased and is no longer the median but is still robust with respect to data outliers.
Lee Pharm owns a bakery and has to decide each morning how many loaves of bread to bake for the day. Any loaves unsold at the end of the day are donated to the local homeless shelter. He would like some help with how to best decide the number of loaves to bake each day. He has available the number of loaves sold for each of the last 13 days:
216, 214, 161, 146, 216, 159, 216, 216, 139, 216, 152, 98, and 116.
In looking at the data, the first thing noticeable is that the number 216 appears quite often. After talking with Lee, it turns out that he typically makes 216 loaves per day, which corresponds to 18 trays of a dozen loaves. Since there was likely unmet demand any day that 216 loaves were sold, after checking with Lee, it turns out that he also has been recording the number of customers asking for bread each day after it was sold out. After adding in this unmet demand, the adjusted potential demand is:
237, 214, 161, 146, 331, 159, 423, 332, 139, 327, 152, 98, and 116.
After confirming with Lee that these past 13 days of demand are reasonably representative of typical demand, the only other information needed for you to provide some guidance is the average selling price of each loaf and the cost to bake each loaf, which turns out to be \$5 and \$1, respectively. Also, Lee agrees that it is reasonable to assume that any demand for a day that is not filled is lost; this is because most customers consume their bread the day of purchase since it is made without preservatives.
A line plot is used because the 13 days of demand represents a time series.
d = [237, 214, 161, 146, 331, 159, 423, 332, 139, 327, 152, 98, 116]
plot(d, legend=false, xlabel="Day", ylabel="Demand", xticks=1:length(d))
Don't see any trend or seasonality that would need to be accounted for in any model of the data. Most importantly, the variability does not show a long-term pattern and so will assume the data is stationary. There are statistical tests for startionarity, but these would not typically be used for so few data points. Since stationarity implies the lack of trend/seasonality, the data can be randomly re-arranged (or permuted), and it should look about the same. We can do this to verify stationarity:
using Random
Random.seed!(2523)
plot(d[randperm(length(d))],legend=false,xlabel="Day",ylabel="Demand",xticks=1:length(d))
plot(d[randperm(length(d))],legend=false,xlabel="Day",ylabel="Demand",xticks=1:length(d))
The two plots of a random permutation of the demands also look stationary, lacking trend/seasonality.
The following determines the daily operating profit if $q$ loaves are baked that morning (i.e., the production rate) and demand turns out to be $d_i$ that day:
$ \quad \begin{eqnarray*} \mbox{Profit: } \pi(q,d_i) &=& \begin{cases} pd_i - cq, & \mbox{if } q > d_i \\ pq - cq = (p - c)\,q, & \mbox{otherwise } \end{cases} \\ &=& p \min \bigl\{ q,d_i \bigr\} - cq \end{eqnarray*}$
p = 5 # Unit price
c = 1 # Unit cost
fπ(q,di) = p*min(q,di) - c*q # Profit per day (fπ ⟹ f\pi<TAB>)
fΠ(q,d) = sum([fπ(q,di) for di in d]) # Total profit (fΠ ⟹ f\Pi<TAB>)
plot(q -> fΠ(q,d), 0, 500, legend=false, xlabel="Production Rate", ylabel="Total Profit")
In order to be able to compare the effectiveness of any procedure, a baseline and upper bound on total profit can be determined. The baseline is the total profit associated with using the current practice of baking 216 loaves per day. An upper bound can be determined by assuming that the number of loaves baked each day exactly matches that day's demand (perfect prediction).
Π_BL = fΠ(216,d)
Π_UB = (p - c)*sum(d)
Π_UB = sum([fπ(di,di) for di in d])
Π_BL, Π_UB, Π_BL/Π_UB
(8517, 11340, 0.7510582010582011)
Since the demand data is assumed stationary, will not try to predict each day's demand; instead, we will assume historical demand data is representative of the variability of any day's demand and will, as a result, use all of the data to determine the same (optimal) number of loaves to bake each day. Note: using 500 as UB for search since it is likely that daily demand might exceed the maximum (423) found using just the 13 days of data.
qᵒ = optimize(q -> -fΠ(q,d), 0, 500).minimizer
display(scatter!([qᵒ],[fΠ(qᵒ,d)]))
qᵒ, fΠ(qᵒ,d), fΠ(qᵒ,d)/Π_UB, fΠ(qᵒ,d)/Π_BL
(330.9999976747946, 9406.999995349586, 0.8295414457980235, 1.104496888029774)
The solution recommends that Lee bake 331 loaves each morning, which should result in an over 10% increase in total profits as compared to baking only 216 loaves. The procedure works because the data is stationary, lacking trend/seasonality. If the data has a trend/seasonality, techniques exist that can try to remove the trend/seasonality so that all of the data can be used; otherwise, for non-stationary data, you would use only the most recent data in any estimate. If, for example, sales have significantly increased in the past week and are expected to remain at this higher level for at least the next few days, then the nonlinear optimization could be rerun using just the last week's data.
Mean and Median: By way of comparison in light of the discussion of MSE (L2) vs. MAD (L1) in NL Opt 2, what the mean or median values for demand were used to determine the number of loaves to bake each morning:
using Statistics
mean(d), fΠ(mean(d),d), median(d), fΠ(median(d),d)
(218.07692307692307, 8541.923076923078, 161.0, 7592.0)
After checking with Lee, the fact that the mean value of 218 is close to the 216 loaf value he was using is no accident: it turns out that he averaged an earlier set of historical demand values to determine the 216 value.
Lee is quite pleased with the prospect of a 10% increase in profits and happens to mention that a large number of his customers preorder their purchases and that the number of orders that he has received each morning seems positively correlated with the number of sales he sees that day. He provides the number of orders available each morning for the past 13 days:
121, 143, 63, 80, 198, 52, 160, 162, 85, 106, 129, 141, and 106During each day, additional walk-in demand occurs, and some orders are canceled; as a result, the final demand at the end of the day differs from the amount ordered at the start of the day. Nevertheless, it is thought that using each day's preorders to help determine the number of loaves to bake that day can increase total profits.
using Plots
d = [237, 214, 161, 146, 331, 159, 423, 332, 139, 327, 152, 98, 116]
o = [121, 143, 63, 80, 198, 52, 160, 162, 85, 106, 129, 141, 106]
display(plot([o d], xlabel="Week", ylabel="Units", xticks=1:length(d),
label=["Orders" "Demand"]))
scatter(o, d, legend=false, xlabel="Orders", ylabel="Demand", label="(o, d)")
The correlation of the two time series can be calculated to see how likely the order data may help in estimating demand:
using Statistics
cor(o,d)
0.5994685344117657
An almost 60% correlation is promising, and so it is likely to be beneficial to include the order data in the analysis. This can be done through linear regression by making the demand a function of order data. Previously, profit was maximized using just the prior demand data. One can think of the demand estimates used previously as single parameter estimation models; linear regression can be used to provide a two-parameter estimation model.
Previousely:
fπ(q,di) = p*min(q,di) - c*q # Profit per day
fΠ(q,d) = sum([fπ(q,di) for di in d]); # Total profit
The total profit function needs to change since a different size ($q_i$) is associated with each different demand ($d_i$). In the single-parameter models considered previously, the same size ($q$) was associated with each different demand ($d_i$).
fΠ(q,d) = sum([fπ(qi,di) for (qi,di) in zip(q,d)]); # Total profit (using ZIP function)
Profit using just order data: What if just the number of orders received by each morning were used to determine the number of loaves to bake that day:
fΠ(o,d)
5969
Linear regression using orders and profit maximization as objective:
fd̂(α,o) = α[1] .+ α[2]o
αΠᵒ = optimize(α -> -fΠ(fd̂(α,o),d), [0. 1.]).minimizer
display(plot!(o -> fd̂(αΠᵒ,o),0,maximum(o), label="Linear"))
fΠ(fd̂(αΠᵒ,o),d), fΠ(fd̂(αΠᵒ,o),d)/Π_UB, fΠ(fd̂(αΠᵒ,o),d)/Π_BL
(9853.636363517973, 0.8689273689169289, 1.1569374619605464)
Question: How would the number of loaves to bake for day 14 be determined if 120 loaves have been ordered by the morning of day 14:
o14 = 120
q14 = Int(round(fd̂(αΠᵒ,o14)))
266
One can consider using other models instead of a line; for example, a quadratic parabola (it is still linear in the coefficients $\alpha$).
fd̂Q(α,o) = α[1] .+ α[2]o .+ α[3]o.^2 # Quadratic model
αΠQᵒ = optimize(α -> -fΠ(fd̂Q(α,o),d), [0. 1. 0.]).minimizer
display(plot!(o -> fd̂Q(αΠQᵒ,o),0,maximum(o), label="Quadratic", legend=:topleft))
αΠQᵒ, fΠ(fd̂Q(αΠQᵒ,o),d)
([-0.07401761718654201 3.5362280660091927 -0.009175237352469189], 9856.37058432405)
The small increase in total profit associated with the quadratic model indicates that it would provide little benefit. Although the 13 data points are too few, with more data, one can perform model validation by randomly selecting a subset of the data to use to train the model (i.e., estimate the coefficients of the model), and then use the remaining data to test what would be resulting total profit using the trained model. In this way, the performance of the linear and quadratic models could be compared.
After further discussion, Lee mentions that he also can access each morning the number of unique visitors to his bakery's website over the past day. We can see if this additinal data can help:
v = [69, 72, 56, 72, 82, 80, 72, 64, 50, 65, 49, 51, 66] # Unique visitors
display(scatter(v, d, legend=false, xlabel="Visitors", ylabel="Demand", label=false))
cor(v,d)
0.46082474901378706
A 46% correlation indicates a moderate relationship that may help. In order to effectively use the visitor data, it can be added to the model along with the order data in a multiple regression:
fd̂M(α,o,v) = α[1] .+ α[2]o .+ α[3]v
αΠMᵒ = optimize(α -> -fΠ(fd̂M(α,o,v),d), [0. 1. 1.]).minimizer
fΠ(fd̂M(αΠMᵒ,o,v),d)
9902.722869055702
Question: How would the number of loaves to bake for day 14 be determined if 120 loaves have been ordered and there have been 45 unique visitors to the website by the morning of day 14:
o14,v14 = 120,45
q14 = Int(round(fd̂M(αΠMᵒ,o14,v14)))
243
| Procedure | Total Profit |
|---|---|
| Baseline | 8517 |
| Single parameter mean (demand) | 8542 |
| Single parameter median (demand) | 7592 |
| Single parameter max profit (demand) | 9407 |
| Using just order data (orders) | 5969 |
| Linear regression (demand; orders) | 9854 |
| Quadratic regression (demand; orders) | 9856 |
| Multiple regression (demand; orders, website visitors) | 9903 |
| Prefect prediction (UB) | 11,340 |
The following dataset contains seven weeks of demand and order data for Lee's bakery:
using DataFrames, CSV, Plots
df = DataFrame(CSV.File("NL_Opt-3-Data.csv"))
first(df,5) # Display only first 5 rows
5 rows × 2 columns
| Demand | Orders | |
|---|---|---|
| Int64 | Int64 | |
| 1 | 434 | 399 |
| 2 | 224 | 128 |
| 3 | 220 | 84 |
| 4 | 361 | 239 |
| 5 | 72 | 22 |
d,o = Int.(df.Demand),Int.(df.Orders)
plot([o d], xlabel="Day", ylabel="Units", label=["Orders" "Demand"])
There are twenty-five switches located throughout an industrial site, and a single monitoring device was just ordered, and it will be located so that it can receive data from the switches. There will be an expensive, dedicated communications line from the device to each of the switches it is monitoring. The coordinates of the switches in the site can be generated by running the following code:
using Random
Random.seed!(73248)
P = 100*rand(25,2)Determine where to locate the device to minimize the length of communications line needed to reach the switches, assuming it can be located anywhere.
using Random, Optim, Statistics, Plots
Random.seed!(73248)
P = 100*rand(25,2)
x0 = mean(P, dims=1)[:] # Using mean switch location as starting point
2-element Vector{Float64}:
55.26100634866732
52.620039507192885
dist2(x,P) = sqrt.(sum((x' .- P).^2, dims = 2))
xᵒ = optimize(x -> sum(dist2(x, P)), x0).minimizer
2-element Vector{Float64}:
54.9883234830379
57.03737843199485
xrng = -10:.1:110
yrng = -10:.1:110
contour(xrng, yrng, (x1,x2) -> sum(dist2([x1, x2],P)))
scatter!(P[:,1], P[:,2], label="Switches")
scatter!([xᵒ[1]], [xᵒ[2]], label="Monitoring Device")
The device has arrived and is ready to be installed, but it turns out that it only has five ports (next time, the specs will be checked before ordering). Luckily, the communications line has not yet been ordered, and so now the problem is to determine where to locate the device to minimize the length of the communications line needed to reach five switches.
function dist2min5pts(x, P)
d = dist2(x, P)
idx = sortperm(d[:])
return d[idx[1:5]]
end
dist2min5pts (generic function with 1 method)
TDᵒ = Inf
iᵒ, xᵒ = NaN, NaN
for i = 1:size(P,1) # Use each point in P as starting point and keep best found
x = optimize(x -> sum(dist2min5pts(x, P)), P[i,:]).minimizer
TD = sum(dist2min5pts(x, P))
println(i,": ",TD)
if TD < TDᵒ
iᵒ,xᵒ,TDᵒ = i,x,TD
end
end
1: 40.66367923059764 2: 82.83992003863787 3: 40.66367922842201 4: 86.40365715634206 5: 40.66367922210874 6: 82.83992004519129 7: 82.8399200393833 8: 40.66367922950188 9: 77.29338760387373 10: 55.28356705439647 11: 40.663679260858395 12: 40.66367922175199 13: 55.28356704909709 14: 81.168947572304 15: 55.28356703113589 16: 40.663679222740186 17: 40.663679223045136 18: 82.83992003666592 19: 55.28356706161319 20: 77.2933876198697 21: 77.29338761229693 22: 40.66367922783404 23: 40.66367922697608 24: 55.28356704383887 25: 40.6636792274368
contour(xrng, yrng, (x1,x2) -> sum(dist2min5pts([x1, x2],P)))
scatter!(P[:,1], P[:,2], label="Switches",
series_annotations = text.([string(i)*" " for i in 1:size(P,1)], :right, 7))
scatter!([xᵒ[1]], [xᵒ[2]], label="Monitoring Device")
idx = sortperm(dist2(xᵒ, P)[:])[1:5] # Five closest switches to device
5-element Vector{Int64}:
25
17
16
3
5
optimize(x -> sum(dist2min5pts(x, P)), P[11,:]) # Nelder-Mead used by default
* Status: success
* Candidate solution
Final objective value: 4.066368e+01
* Found with
Algorithm: Nelder-Mead
* Convergence measures
√(Σ(yᵢ-ȳ)²)/n ≤ 1.0e-08
* Work counters
Seconds run: 0 (vs limit Inf)
Iterations: 43
f(x) calls: 84
optimize(x -> sum(dist2min5pts(x, P)), P[11,:], LBFGS()) # Compare L-BFGS
* Status: success
* Candidate solution
Final objective value: 4.066368e+01
* Found with
Algorithm: L-BFGS
* Convergence measures
|x - x'| = 5.56e-06 ≰ 0.0e+00
|x - x'|/|x'| = 6.02e-08 ≰ 0.0e+00
|f(x) - f(x')| = 9.50e-11 ≰ 0.0e+00
|f(x) - f(x')|/|f(x')| = 2.34e-12 ≰ 0.0e+00
|g(x)| = 4.57e-09 ≤ 1.0e-08
* Work counters
Seconds run: 1 (vs limit Inf)
Iterations: 13
f(x) calls: 32
∇f(x) calls: 32
Warning: Automatic differentiation does not work for this problem. It could be because some of the objective function code is not pure Julia at a low level (?), but beware that the solution is reported as being a "success"!
optimize(x -> sum(dist2min5pts(x, P)), P[11,:], NewtonTrustRegion(), autodiff=:forward)
* Status: success
* Candidate solution
Final objective value: 9.983148e+01
* Found with
Algorithm: Newton's Method (Trust Region)
* Convergence measures
|x - x'| = 0.00e+00 ≤ 0.0e+00
|x - x'|/|x'| = 0.00e+00 ≤ 0.0e+00
|f(x) - f(x')| = NaN ≰ 0.0e+00
|f(x) - f(x')|/|f(x')| = NaN ≰ 0.0e+00
|g(x)| = NaN ≰ 1.0e-08
* Work counters
Seconds run: 0 (vs limit Inf)
Iterations: 0
f(x) calls: 1
∇f(x) calls: 1
∇²f(x) calls: 1