4.3: Production-Inventory Planning: Single Product
ISE 754: Logistics Engineering, Fall 2026
The interdependence of production decisions quickly overwhelms simple planning rules. For these complex problems, reliable heuristics are rarely available. MILP provides the essential framework to systematically evaluate every trade-off and identify the optimal plan.
No new Julia packages used.
No new Logjam functions used.
Companion script
4-net-3.jl, runnable Julia extracted from the lecture’s code blocks.
1. Production-inventory planning
Production-inventory planning models are a primary application of mathematical programming in industry. They provide a means to make complex decisions over a rolling planning horizon. Decisions are complex because each has an impact across multiple time periods and stages in the production process. Reliable heuristics are rarely available. A time-expanded network is used to model multi-period production-inventory problems. Models are re-solved each time period, using the latest mix of firm and forecasted orders over the planning horizon.
Among the Networks models, production-inventory planning stands somewhat apart, because it leaves the geographic representation of a problem behind. A production-inventory problem has two dimensions: time, since production is planned over several periods, and the stages of production, the output of each stage feeding the next until the last sends the product to the customer. The stages need not all occur in the same period, and when they do not, the product is carried as inventory at a real cost, to be weighed against the cost of production. What makes the problem interesting is that the resources are capacitated, which makes it a tactical problem rather than a strategic one: how much capacity to build is not being determined, and the problem is how best to deploy the capacity available to produce a series of products. Once setup costs (Sec. 3) or several products (Lecture 4.4) are part of it, it is the one point in the course where a mixed-integer linear program is really needed, since there are no good heuristics for the problem and no other way of looking at it. The uncapacitated facility location (UFL) problem was also solved as a mixed-integer linear program (Lecture 2.7), but there it was more than was needed, since a good heuristic exists for it (Lecture 2.4), and a mixed-integer linear program is worth using for the UFL only when constraints are added that a heuristic would find hard to incorporate.
Production-inventory planning is the core of the course: if the course covered nothing else, this would be the one thing to cover.
NoteHow planning software grows
Much of it is embedded in the software used in logistics engineering, so only a developer at one of the companies that build the software works with the models directly. A field engineer, dealing with a client’s problems, still has to know enough of it to judge whether a problem, a constraint, or some other aspect of the client’s supply chain or logistics network can be modeled. Sometimes the model already has the capability and it is easily done. The interesting cases are those where it does not, and the client is large enough that the capability is developed for it. That is how such software evolves: nobody designs it from scratch; it grows one project after another, each slightly different from the last, and a solution to one client’s particular problem is easy to generalize for the next. After ten to twenty years of this, a commercial product has its full feature set.
A single product, with a model that handles that product alone, is the place to start; Lecture 4.4 treats multiple products that share some of the resources, where how they share them is determined within the model.
Decisions in these models are complex because each one affects several periods and stages. A production-inventory model is needed only when a decision made now limits what can be done later: with weekly periods, a decision made this week can rule out something three weeks from now. A system in which everything happens within a couple of hours, and in which any decision can be redone, never needs one. The models suit production that is committed once it is set up, such as a refinery set up for one kind of product, and resources with long scheduled downtime: a two-month maintenance shutdown is planned years ahead, and since the low-demand season is too short to maintain every resource in it, when each one is maintained is itself a decision. Decisions like these soon outrun any feel for a good answer, which leaves a planner more reliant on the models than in most problems. Road networks and location problems can be seen on a map; production-inventory problems cannot be visualized well, and the mathematical program is the best way of looking at them.
Demand is handled differently here than in the rest of the course. The earlier lectures assumed a constant demand, as strategic problems typically do, since they seldom need to incorporate variability at all, and two of the Inventory lectures (Lectures 6.1 and 6.2) model variability explicitly. Production-inventory planning is a middle ground: its models are deterministic, but they handle variability by having a different demand in each period over the planning horizon. That handles both the variability and the need for a deterministic model, without stochastic programming or simulation, neither of which can handle models as complex as production-inventory models can be.
Lecture 1.3 also assumed a constant average demand, and captured variability through the factory-physics equations, among them the VUT (variability, utilization, time) equation for cycle time (Eq. 6 in Lecture 1.3). That is a different type of variability from the one considered here. The variability here is structural: it is deterministic in each period, and the different demand in each period captures it.
Although technically the capacity is variable, it is really not variable: it is either available or not available. In Sec. 2, the decision about whether capacity is available or not is made outside of the model. Prior to effective mixed-integer linear programming, this was the only way these kinds of problems could be solved. A logistics engineer would have to make the capacity decisions, and capacity that is not available would just be set to zero. Then this would be solved as a linear program, which was doable even in the 1980s. Making the decision internally to the model, the subject of Sec. 3 on setup costs, was only possible when mixed-integer linear programming capability came online in the early to mid 1990s. The setup costs were also there in the problems of Sec. 2, but the cost of the setups was not needed because the setups had already been made outside the model. This is an important distinction, and the whole point of this lecture.
Neither the planning horizon nor the demand over it is given to a logistics engineer. As with many things, determining an appropriate planning horizon is the engineer’s responsibility, and part of the whole process of engineering a solution to a production-inventory problem. There is not necessarily a single answer: different engineers can arrive at different horizons, some better than others. Forecasting the demand over the horizon is part of the same responsibility. The examples of Secs. 2 and 3 are given their horizon and their demand, so that the focus is on the models; Sec. 4 then takes up where the demand comes from, and Sec. 5 how long the horizon should be.
The time-expanded network of Fig. 1 is the basic setup of every production-inventory model, and understanding it is most of what is needed to understand the models. It is the network of a single product. Each row of nodes is a stage of production, two in the figure, and each column a period being planned, typically a week or a month; a year would usually be too long. The most important input is the demand forecast for each period. The current period’s demand, 20 tons in the figure, is probably all orders already received, so it is a very good number; the second period’s is perhaps two-thirds booked orders and one-third forecast, the third’s one-third booked and two-thirds forecast, and further out more and more of it is forecast, until all of it is. With monthly periods the model is re-solved every month with the actual orders and the best forecast of the rest, and each solution says how production is to run, above all in the current period: how much to produce at each stage, how much of it meets demand, and how much is carried into the next period as inventory. Everything is measured in tons. The x variables are the tons produced, and the y variables the tons of inventory, the tons at a stage carried from one period into the next, where they can go on to the next stage and ultimately meet demand, or be carried again.
Show the code that draws this figure
# Code block 1: drawing the time-expanded network of a production plan# Stage i is row i and period j column j; each arc runs from node outline# to node outline, and each label is a symbol in black with its value in# red, as on the slide the figure is redrawn from; given a plan's# production X and inventory Y, their values replace the capacities and# the inventory symbols, in blueink, datac = colorant"#252525", colorant"#d21f26"flowc = colorant"#1f5fa8"# a plan's flows# Every symbol in the font of the lecture's equations, New Computer# Modern: italic for the letter, upright for its digits and operatorsmi = Makie.MathTeXEngine.texfont(:italic)mr = Makie.MathTeXEngine.texfont(:math)it(s) =rich(s; font = mi) # a symbol, in italicsup(s; kw...) =rich(s; font = mr, kw...) # digits, operatorsmth(s) =rich([isletter(c) ? it(string(c)) :up(string(c)) for c in s]...) # a mixed termsym(b, sub, sup =nothing) = sup ===nothing ? # b with subscript subrich(it(b), subscript(up(sub))) :# and superscript suprich(it(b), subsup(up(sub), mth(sup)))num(v) = (r =round(v, digits =2); isinteger(r) ? Int(r) : r)val(v; c = datac) =up(" = $(num(v))"; color = c) # a valuefunctiondrawpinet(cᵖ, u, cⁱ, d, y⁰, yᵗ⁺¹; X =nothing, Y =nothing, unit ="month", w =680, h =500, r =0.17, ls =16, dx =2.3, dy =2.0, g =0.07) # g: label to its arc m, t =length(cᵖ), length(d) fig =Figure(size = (w, h)) ax =Axis(fig[1, 1]; aspect =DataAspect())hidedecorations!(ax); hidespines!(ax)P(i, j) = ((j -1) * dx, -(i -1) * dy) # stage i, period jarc(p, q) =arrows2d!(ax, [p[1]], [p[2]], [q[1] - p[1]], [q[2] - p[2]]; color = ink, shaftwidth =1.3, tipwidth =8, tiplength =8)lbl(p, s; al = (:left, :center), rot =0.0, off = (0, 0)) =text!(ax, p[1], p[2]; text = s, color = ink, fontsize = ls, align = al, rotation = rot, offset = off) ℓ, ℓs =1.0, 1.5# end arcs, supply arcsfor j in1:t x0, y0 =P(1, j) # supply into stage 1arc((x0, y0 + ℓs), (x0, y0 + r))for i in1:m # production, stage i xa, ya = i ==1 ? (x0, y0 + ℓs) :P(i -1, j) .- (0, r) ym = (ya +P(i, j)[2] + r) /2lbl((xa +0.1, ym +0.2),rich(sym("c", "$i", "p"), val(cᵖ[i])))lbl((xa +0.1, ym -0.2), X ===nothing ?rich(sym("x", "$i,$j"), up(" ≤ "), sym("u", "$i,$j"),val(u[i, j])) :rich(sym("x", "$i,$j"), val(X[i, j]; c = flowc))) i >1&&arc((xa, ya), P(i, j) .+ (0, r))end xm, ym =P(m, j) # demand out of stage marc((xm, ym - r), (xm, ym - ℓ))lbl((xm, ym - ℓ -0.22), rich(sym("d", "$j"), val(d[j])); al = (:center, :center))endfor i in1:m # inventory arcs xa, ya =P(i, 1) .- (ℓ +0.4, 0) # initial xb =P(i, 1)[1] - rarc((xa, ya), (xb, ya))lbl(((xa + xb) /2, ya + g),rich(sym("y", "$i,1"), up(" = "), sym("y", "$i", "0")); al = (:center, :bottom))lbl(((xa + xb) /2, ya - g), val(y⁰[i]); al = (:center, :top))for j in1:t-1 p, q =P(i, j) .+ (r, 0), P(i, j +1) .- (r, 0)arc(p, q)lbl(((p[1] + q[1]) /2, p[2] + g), # cost, just aboverich(sym("c", "$i", "i"), val(cⁱ[i])); al = (:center, :bottom))# flow, below; a y alone has no tall glyph, so it is raised by# about half its font size, which puts its top as close to its# arc as the tops of the values below the other arcslbl(((p[1] + q[1]) /2, p[2] - g), Y ===nothing ?sym("y", "$i,$(j+1)") :rich(sym("y", "$i,$(j+1)"), val(Y[i, j+1]; c = flowc)); al = (:center, :top), off = Y ===nothing ? (0, 0.52ls) : (0, 0))end p =P(i, t) .+ (r, 0) # final q = p .+ (ℓ +0.4- r, 0)arc(p, q)lbl(((p[1] + q[1]) /2, p[2] + g),rich(sym("y", "$i,$(t+1)"), up(" = "), sym("y", "$i", "t+1")); al = (:center, :bottom))lbl(((p[1] + q[1]) /2, p[2] - g), val(yᵗ⁺¹[i]); al = (:center, :top))endfor i in1:m, j in1:t # nodes, over the arcspoly!(ax, Circle(Point2f(P(i, j)...), r); color =:white, strokecolor = ink, strokewidth =1.8)end xc = (t -1) * dx /2# what the arcs arelbl((xc, ℓs +0.15), "Supply"; al = (:center, :bottom)) yb =-(m -1) * dy - ℓ -0.55lbl((xc, yb), "Demand"; al = (:center, :top))for j in1:tlbl((P(1, j)[1], yb -0.4), "$j"; al = (:center, :top))endlbl((xc, yb -0.75), "Period ($(unit))"; al = (:center, :top)) xl =-(ℓ +0.4) -0.25lbl((xl, -(m -1) * dy /2), "Initial inventory"; al = (:center, :bottom), rot =π/2) xr = (t -1) * dx + ℓ +0.4+0.25lbl((xr, -(m -1) * dy /2), "Final inventory"; al = (:center, :top), rot =π/2)for i in1:mlbl((xl -0.55, P(i, 1)[2]), "$i"; al = (:center, :center))endlbl((xl -0.85, -(m -1) * dy /2), "Stage"; al = (:center, :bottom), rot =π/2)return figenddrawpinet([200, 800], [50050; 6000], [5, 25], [20, 10, 15], [0, 0], [35, 0])
Figure 1: The time-expanded network of a single product made in two stages over three monthly periods, with the data of Ex. 3. Production enters each stage’s node at a cost of c_i^p per ton, up to the stage’s capacity u_{ij} in the period; inventory y_{ij} is carried from period j-1 to period j at a cost of c_i^i per ton; and the last stage meets each period’s demand d_j. Everything in black is a symbol, and everything in red its particular value. The horizontal arcs, left to right, are the inventory flows, and the vertical arcs, top to bottom, are the production flows.
2. Capacity determined outside the model
For a single product without setups, a plan can usually be found by hand, without a linear program: every arc of the time-expanded network points forward, so each ton can be made as late as capacity allows, working back from the last stage, by the procedure Model 2 states below. The procedure is exact only under its assumptions. A storage limit, an upper bound on the amount of inventory that can be carried from one period to the next, can break it: in Ex. 4 the latest plan would hold more at a stage than it can store, and only the linear program finds the plan that does not. Setup costs (Sec. 3) make which periods to produce in a decision of its own, which takes a mixed-integer linear program. The linear program, Model 1, is the model that carries these; the procedure is how a plan is found, and checked, by hand.
With several products (Lecture 4.4), the problem is formulated as an integer program, and there is no simple heuristic to rely on; that is where the solver is really needed.
For production-inventory problems in general, no heuristic is known short of solving the whole problem, because a decision made now affects what happens not only now but four, five, six, seven or eight periods later.
2.1 Linear program
The following production-inventory LP planning model is for a single product, produced in m stages with a rolling horizon of t periods (e.g., months). Determining if production capacity u_{ij} should be available for each period, e.g., u_{ij} = u_i or u_{ij} = 0 is used to plan production.
The first model is for a production system that produces a single product, and it makes no internal decision about what to produce when. What makes the problem more than trivial is that the resources are typically shared with other products, and the model stays a linear program by handling that through the capacities u_{ij}, the production capacity of stage i in period j. In a period in which a resource is available for the product, all of its capacity is, say 10 tons a month; in other periods it makes other products. Which product to make is not determined by the model, since that is the multiple-product problem of Lecture 4.4; instead, the person operating the model sets the capacity to zero in each period in which the resource is not available. That is very labor-intensive, and with multiple products today’s solvers can simply solve the problem with that decision in it. Twenty or thirty years ago, before solver performance improved enormously, there was no other way than solving a lot of linear programs: the logistics engineer determined when to make capacity available, and the linear program then found, for those decisions, the best flow of product and how much of it to carry over as inventory, which still took a great many of the details off the engineer’s hands. Modeling the decision directly takes it away from the engineer altogether and makes still more complex problems possible. In the model, the objective is the unit production cost and the inventory carrying cost; the variables x_{ij} are how much to produce at each stage in each period, continuous since the model is a linear program; and the variables y_{ij} are the inventory carried from period to period. The superscript p marks a production cost and i an inventory cost. The network of Fig. 1 is the easiest way to understand the formulation.
Fig. 1 shows the model on its data, a network flow problem. Raw materials go into stage 1, and the first column of nodes is the first of three months. The arc into stage 1 in month 1 has an upper bound, its capacity u_{1,1}, here 50 tons, and a cost of $200 per ton; u_{1,2} is stage 1’s capacity in month 2, and so on. The initial inventory, y_1^0 at stage 1, is where the plan stands now, carried over from the previous run, since the linear program is re-solved every month: next month’s run makes month 2 the first month and extends the plan one month further, and the inventory on the horizontal arcs becomes the initial inventory of the new first month. Because the model is not run to infinity, the final inventory is the inventory wanted when the plan ends, and typically it is not zero: a final inventory of zero would not be realistic for an ongoing production process that is not being shut down in three months. An ongoing operation would also typically have a positive initial inventory; Fig. 1 sets it to zero for simplicity. At stage 2, production costs $800 per ton, with a capacity of 60 tons, and the flow out of stage 2 is the demand, 20, 10 and 15 tons in the three months. Carrying a ton from month 1 to month 2 costs $25 at stage 2, more than at stage 1, because $800 more has been spent processing it. Some capacities are zero, stage 1’s in month 2 and stage 2’s in months 2 and 3, because the resource may be making another product then, or be down for scheduled maintenance.
What the stages mean in context is the steps of production, since production does not all happen at one point. Stage 1 might take the raw material, a metal or a plastic, and cast it into a particular form, and stage 2 take the casting and process it further, say trimming features off its edges that need correcting. A stage can also be the physical transportation of the product: the processing of one stage takes place at one location, and the product is then put on a plane and flown to another country for the next stage, perhaps a second operation, and then flown somewhere else again, so that a complex production process can easily have dozens of stages. The inventory cost reflects this. The product at stage 2 is in a sense the same product as at stage 1, but not really: $200 has been spent to produce a ton of it at stage 1, and another $800 at stage 2. Carrying a ton costs the carrying rate times what has been invested in it, so a ton at stage 2, carrying what was spent at both stages, costs more to carry than a ton at stage 1; Ex. 3 calculates the two costs, the $5 and $25 per ton of Fig. 1, in Eq. 2. In essence, each stage adds value. That is obvious in manufacturing, which changes the physical form of a product, but it applies as well to a product put in a shipping container in China when $3,000 is spent to transport the container to the port of Long Beach: the $3,000 is added to the value of all the products in the container.
So the product is worth more when it is taken off the ship in Los Angeles than when it left China. The difference is often small enough to ignore, but every manufacturing or transportation step adds value, even one that leaves the product physically unchanged, and that is how the inventory costs of Fig. 1 are to be read. The solver makes no use of the network’s structure; the model is solved simply as a linear program. For a single product, as here, it is a minimum cost network flow problem, Model 3 in Lecture 4.2. Once integer or binary variables are added (Sec. 3), it is no longer the classic minimum cost network flow problem, but as a linear program it is exactly that.
The question is how to do this at least cost. The two costs that can be minimized are the cost of production and the cost of inventory. At each node, the inventory carried in and what is produced are the input, and the output is either sent on to the next stage, used to satisfy demand, or carried over as inventory. Model 1 states the model in words and as a linear program.1
minimize: total production and inventory cost over the planning horizon
solve for:
(a) production at each stage in each period;
(b) inventory at each stage carried into each period.
subject to:
(a) flow balance at each stage that feeds the next: what is produced or carried in is either passed on to the next stage or carried out;
(b) flow balance at the last stage, which feeds demand: what is produced or carried in either meets the period’s demand or is carried out;
(c) capacity: production at each stage in each period is no more than its capacity;
(d) initial and final inventory: inventory at each stage at the start and at the end of the horizon is given;
(e) storage: inventory at each stage between periods is no more than its storage capacity.
assumptions:
(a) production cost and inventory cost per ton at each stage are constant;
(b) demand in each period is given;
(c) production can pass through every stage within the period it starts in.
=inventory at stage i carried from period j-1 to period j (ton)
c_i^p
=production cost at stage i ($/ton)
c_i^i
=inventory cost at stage i ($/ton-period)
=h\,(c_1^p + c_2^p + \cdots + c_i^p)
h
=inventory carrying rate per period ($/$-period)
d_j
=demand in period j (ton)
u_{ij}
=capacity of stage i in period j (ton)
u_i^y
=storage capacity at stage i (ton)
y_i^0
=initial inventory at stage i (ton)
y_i^{t+1}
=final inventory at stage i (ton)
m
=number of production stages
t
=number of periods in the planning horizon.
Both lines labeled (d) in Eq. 1 are row (d) on purpose: one constraint, the inventory given at the start and at the end of the horizon.
In the objective function, inventory costs are summed from 2 to t + 1 because the period 1 initial inventory cost was accounted for in period 0 prior to period 1 (rolling horizon). Constraints (a) and (b) in the model are the flow-balance equations for each node, while the remaining constraints represent lower- and upper-bounds on the decision variables x and y.
The inventory sum in the objective runs from 2 to t + 1 for a reason. There are t + 1 inventory variables at each stage, four in the three periods of Fig. 1, but summing their costs from 1 to t + 1 would double count. In a rolling horizon, this month’s plan carries those four, and next month’s cuts off the first period, adds one at the end, and solves the same model again. With the cost counted every month, the inventory carried from one plan into the next would appear in both, as the final inventory of one and the initial inventory of the other, and over time the plans would show an inventory cost that is not actually incurred. Leaving out the cost of the initial inventory avoids that double counting. The final inventory could have been left out instead, and the initial included: one or the other, but not both. That is an important aspect of these models, and it is why they are written with indices rather than set notation. After the objective function, the model is just flow balance: at every stage but the last, everything sums to zero. The only exception is the last stage, whose sum equals the negative of the demand, because in the net supply formulation of the minimum cost network flow problem a demand is a negative net supply. So Model 1 is another example of the minimum cost network flow problem of Lecture 4.2, a specialization of it.
Because of the need to sum over different portions of the data, indices are used in place of set notation in production-inventory planning models; for example, while j \in T could be used in place of j = 1, \dots, t, it would be difficult to use set notation to represent i = 1, \dots, m-1.
The flow-balance constraints are the only constraints in the model because the capacity constraints are handled as upper bounds on x_{ij}.
Model 1 implementation: Production-inventory planning
# Model: production-inventory planningusingJuMP, HiGHS, Logjamfunctionprodinv(cᵖ, cⁱ, d, u, LBʸ, UBʸ) m, t =length(cᵖ), length(d) model =Model(HiGHS.Optimizer)@variable(model, 0≤ x[i=1:m, j=1:t] ≤ u[i,j]) # (c)@variable(model, LBʸ[i,j] ≤ y[i=1:m, j=1:t+1] ≤ UBʸ[i,j])@objective(model, Min,sum(cᵖ[i]*x[i,j] for i=1:m, j=1:t) +sum(cⁱ[i]*y[i,j] for i=1:m, j=2:t+1))@constraint(model, flowbal[i=1:m-1, j=1:t], # (a)-x[i,j] + x[i+1,j] - y[i,j] + y[i,j+1] ==0)@constraint(model, laststage[j=1:t], # (b)-x[m,j] - y[m,j] + y[m,j+1] ==-d[j])set_silent(model)optimize!(model)iftermination_status(model) == JuMP.OPTIMALreturn (TC =objective_value(model), x =snapvals(value.(x)), y =snapvals(value.(y)))enderror(termination_status(model))end
1
Rows (d) and (e) are bounds on y: its lower bound LBʸ and upper bound UBʸ are the same initial inventory in period 1 and the same final inventory in period t + 1, which fixes both, and between them zero and the storage capacity.
prodinv (generic function with 1 method)
2.2 Flows determined by flow balance
Flow balance makes it possible to determine some of the flows in a production-inventory network by hand. Given a diagram of the network with only some of its flows, the missing flows are filled in from the flow balance at its nodes. Doing so is a good way of demonstrating an understanding of the basic structure of these networks.
Missing flows are not something to expect in practice, and filling them in is not part of the normal use of these models. Where some are missing, this is a way to fill them in; it works only because the structure of the model is such that some of the data can be determined from the rest. It is brought up here because it is a useful way of testing an understanding of that structure. Were flows to be found to any large extent, it would be coded up, so that the model determined them automatically. Exs. 1 and 2 are two such tests, the first with one flow missing and the second with data of several kinds missing.
Example 1: One missing flow
Determine the missing value of flow c in the network of Fig. 2(a).
Show the code that draws this figure
# Code block 2: Ex. 1's network as given, and its missing flow found# frame! draws the flows N on the network, each marked by what it is in# N₀, the network as given: given red, unknown black and named, found blue# and named; a node in hl is circled in blue. A flow is named by its array# and indices: (:X, i, j) production, (:Y, i, j) inventory, (:d, j) demandfunctionframe!(ax, N, N₀, name; hl = [], r =0.16, ls =12) m, t =size(N.X)hidedecorations!(ax); hidespines!(ax)P(i, j) = ((j -1) *1.7, -(i -1) *1.6) # stage i, period jarc(p, q) =arrows2d!(ax, [p[1]], [p[2]], [q[1] - p[1]], [q[2] - p[2]]; color = ink, shaftwidth =1.0, tipwidth =6, tiplength =6)functionlab(k) # a flow's label, color v, v₀ = N[k[1]][k[2:end]...], N₀[k[1]][k[2:end]...]ismissing(v₀) ||returnstring(v₀), datac s =rich(it(name[k]), up(ismissing(v) ? " = ?":" = $v"))return s, ismissing(v) ? ink : flowcendtx(p, k, al) = ((s, c) =lab(k); text!(ax, p[1], p[2]; text = s, color = c, fontsize = ls, align = al))for j in1:t # production and demand x0, y0 =P(1, j)arc((x0, y0 +0.9), (x0, y0 + r))tx((x0 +0.08, y0 +0.5), (:X, 1, j), (:left, :center))for i in2:m p, q =P(i -1, j), P(i, j)arc(p .- (0, r), q .+ (0, r))tx((p[1] +0.08, (p[2] + q[2]) /2), (:X, i, j), (:left, :center))end xm, ym =P(m, j)arc((xm, ym - r), (xm, ym -0.7))tx((xm, ym -0.8), (:d, j), (:center, :top))endfor i in1:m # inventory p =P(i, 1) .- (0.75, 0)arc(p, P(i, 1) .- (r, 0))tx(p .- (0.08, 0), (:Y, i, 1), (:right, :center))for j in1:t-1 p, q =P(i, j) .+ (r, 0), P(i, j +1) .- (r, 0)arc(p, q)tx(((p[1] + q[1]) /2, p[2] +0.12), (:Y, i, j +1), (:center, :bottom))end q =P(i, t) .+ (0.75, 0)arc(P(i, t) .+ (r, 0), q)tx(q .+ (0.08, 0), (:Y, i, t +1), (:left, :center))endfor i in1:m, j in1:t # nodes, over the arcs h = (i, j) in hlpoly!(ax, Circle(Point2f(P(i, j)...), r); color =:white, strokecolor = h ? flowc : ink, strokewidth = h ? 2.6:1.4)endxlims!(ax, -1.45, 5.3); ylims!(ax, -2.85, 1.0)return axendX = [9019; 19012] # production at stage i in period jd = [7, 12, 8] # demand in period jY = [10007; # inventory at stage i into period0missing04] # j, missing where unknownYᵒ =copy(Y)Yᵒ[2, 2] = Y[2, 1] + X[2, 1] - d[1] # flow balance at stage 2, month 1N₀ = (X = X, Y = Y, d = d)N = (X = X, Y = Yᵒ, d = d)fig =Figure(size = (680, 260))for (k, (Nk, hl, title)) inenumerate([(N₀, [], "(a) As given"), (N, [(2, 1)], "(b) Balance at stage 2, month 1")]) ax =Axis(fig[1, k]; aspect =DataAspect(), title = title, titlesize =14, titlefont =:regular)frame!(ax, Nk, N₀, Dict((:Y, 2, 2) =>"c"); hl = hl, ls =15)endcolgap!(fig.layout, 8)fig
Figure 2: Ex. 1’s missing flow found by flow balance: (a) the network as given, its production, demand and inventory flows in red and the unknown inventory flow c; (b) the node of stage 2 in month 1 balanced, circled in blue, and the flow it determines in blue.
At stage 2 in the first month, 19 tons come in, and 7 are needed for the demand, which leaves 12, and the only place the 12 can go is into the second month as inventory, so c is 12; every other node nets the rest to zero. Flow balance always holds, so it is a way of validating, whenever a plan has any complexity, that things are working the way they should: the Balance check of Lecture 1.2.
c = 12.
NoteFind missing flows using Julia
Code block 3 finds missing flows the way they are found by hand, node by node: fillflows balances each node that has exactly one unknown flow, records the step, and repeats until no flow is missing. Production x_{ij} is the matrix X, inventory y_{ij} the matrix Y, and demand d_j the vector d, each with missing for an unknown, and the flows found are returned in N:
# Code block 3: missing flows found node by node by flow balance# At node (i, j), inventory in plus production equals what goes on, to# stage i + 1 or to demand, plus inventory out; a node with one unknown# flow determines it, and each step records the node and the flow foundfunctionfillflows(X, Y, d) m, t =size(X) N = (X =copy(X), Y =copy(Y), d =copy(d))v(k) = N[k[1]][k[2:end]...] steps =Tuple[]whileany(A ->any(ismissing, A), N) n =length(steps)for i in1:m, j in1:t F = [(+1, (:Y, i, j)), (+1, (:X, i, j)), # in (-1, i < m ? (:X, i +1, j) : (:d, j)), # on (-1, (:Y, i, j +1))] # out u = [f for f in F if ismissing(v(f[2]))]length(u) ==1||continue s, k = u[1] N[k[1]][k[2:end]...] =-s *sum(s′ *v(k′) for (s′, k′) in F if k′ != k)push!(steps, (i, j, k))endlength(steps) > n ||error("flow balance cannot find the rest")endreturn N, stepsendX = [9019; 19012] # production at stage i in period jd = [7, 12, 8] # demand in period jY = [10007; # inventory at stage i into period0missing04] # j, missing where unknownN, steps =fillflows(X, Y, d)@show stepsN.Y[2, 2]
1
A flow is named by its array and its indices, (:Y, 2, 2) for y_{2,2},(:X, 1, 3) for x_{1,3} and (:d, 1) for d_1, and v(k) reads the flow named k.
2
The unknown flow’s sign s is +1 if it comes into the node and -1 if it goes out. Flow balance, what comes in less what goes out equal to zero, gives the flow as -s times the signed sum of the node’s other three.
3
n counts the steps before a pass. If a pass finds nothing, the rest cannot be found by flow balance alone, and the error says so.
steps = Tuple[(2, 1, (:Y, 2, 2))]
12
fillflows finds the one step taken by hand, at stage 2 in month 1, and the same c.
The flows behave like water flowing through a pipe, and the analogy is exact. With some of the information missing, the rest can be inferred from these two ideas, which is how the problem is done by hand.
Example 2: Unknowns of several kinds
Determine the missing values p,q,r and s in the network of Fig. 3(a): a production flow, an inventory flow, a final inventory and a demand.
Show the code that draws this figure
# Code block 4: Ex. 2's missing data found by fillflows, a frame a stepX = [40030; missing2022] # production at stage i in period jd = [20, missing, 17] # demand in period jY = [1025missingmissing; # inventory at stage i into period611712] # j, missing where unknownN₀ = (X = X, Y = Y, d = d)N, steps =fillflows(X, Y, d)name =Dict((:X, 2, 1) =>"p", (:Y, 1, 3) =>"q", (:Y, 1, 4) =>"r", (:d, 2) =>"s")done = [s[1:2] for s in steps] # the nodes balanced, in orderleft = [(i, j) for i in1:2 for j in1:3 if (i, j) ∉ done]fig =Figure(size = (680, 680))Nk = (X =copy(X), Y =copy(Y), d =copy(d)) # the data known so farfor k in0:length(steps) +1if k ==0 hl, title = [], "(a) As given"elseif k <=length(steps) i, j, f = steps[k] Nk[f[1]][f[2:end]...] = N[f[1]][f[2:end]...] hl, title = [(i, j)], "($('a'+ k)) Balance at stage $i, month $j"else# the nodes left over check hl, title = left, "($('a'+ k)) Check: the nodes left over"end ax =Axis(fig[k ÷2+1, k %2+1]; aspect =DataAspect(), title = title, titlesize =13, titlefont =:regular)frame!(ax, Nk, N₀, name; hl = hl, ls =14)endcolgap!(fig.layout, 8); rowgap!(fig.layout, 4)fig
Figure 3: Ex. 2’s missing data found by flow balance, one frame per step of fillflows (Code block 3): (a) the network as given, its known flows in red and four unknowns, the production flow p, the inventory flow q, the final inventory r and the demand s; (b) to (e) one node balanced at each step, circled in blue, and the value it determines in blue; (f) the two nodes left over, which balance too.
Stage 1 is balanced first, month by month. In month 1, the 10 tons of initial inventory and the 40 produced come in, and 25 are carried into month 2, so what goes on to stage 2 is p = 10 + 40 - 25 = 25 (frame b). In month 2, the 25 carried in and nothing produced feed the 20 that stage 2 makes, which leaves q = 25 + 0 - 20 = 5 to carry into month 3 (frame c). In month 3, those 5 and the 30 produced feed the 22 that stage 2 makes, which leaves the final inventory, r = 5 + 30 - 22 = 13 (frame d), so r could be found only after q. At stage 2, month 2’s node gives the demand: the 11 carried in and the 20 made, less the 7 carried out, s = 11 + 20 - 7 = 24 (frame e). The two nodes left over check the answer (frame f): in month 1, 6 + 25 comes in and 20 + 11 goes out, and in month 3, 7 + 22 comes in and 17 + 12 goes out. The first of the two is also a second route to p, since its only unknown was p:p = 20 + 11 - 6 = 25, the Triangulate check of Lecture 1.2.
p = 25, q = 5, r = 13, s = 24.
Not all missing variables can be recovered: two unknowns that meet at only one node are known only through their sum, as the demand in month 3 and stage 2’s final inventory would be in Ex. 2, and fillflows then stops with its error.
2.3 Greedy procedure
Inventory costs more to carry at each later stage, because more has been spent producing the product by then, so it often makes sense to carry inventory at an early stage, such as stage 1, and to produce each ton only when it is actually needed.
For a single product, something like a simple greedy procedure can be used, because the network is essentially feed forward, and solving the problem as a linear program can be overkill.
Every arc of the time-expanded network points forward, to a later period or to the next stage, so a plan can be found by working back along it, from the last stage: the last stage’s latest plan leaves every earlier stage the most time to supply it. Model 2 states the procedure, which is how a plan is found by hand. It is a feasibility model, not an optimization model, in the terms of the Model Format Reference: it determines the production with no preference among the plans that work, so it has no objective function, and nothing in it is chosen by cost: under its assumption (a), the production cost is the same for every plan that satisfies the demand, which leaves only the inventory cost, and under the rest, producing each ton as late as capacity allows is what keeps that least. The rule takes the place of the objective, and Ex. 3 checks that its plan costs what the linear program’s does.
solve for:
(a) production at each stage in each period.
subject to:
(a) as late as possible: each stage produces its requirement in the latest periods its capacity allows, the last stage’s requirement being the demand and each other stage’s the production of the stage after it;
(b) initial and final inventory: a stage’s initial inventory meets its earliest requirements, and its final inventory is produced as late as the rest.
assumptions:
(a) production cost and inventory cost per ton at each stage are the same in every period;
(b) inventory cost does not decrease from one stage to the next;
(c) no storage limit binds;
(d) production can pass through every stage within the period it starts in.
Model 2: Greedy as-late-as-possible procedure
Model 2 formulation: Greedy as-late-as-possible procedure
algorithm latestplan;{ input: capacities uᵢⱼ, demands dⱼ, initial and final inventories }{ yᵢ⁰ and yᵢᵗ⁺¹, stages i = 1 … m, periods j = 1 … t }{ output: production xᵢⱼ and inventory yᵢⱼ, or infeasible }begin for i := m downto 1 do begin yᵢ₁ := yᵢ⁰; yᵢ,ₜ₊₁ := yᵢᵗ⁺¹; for j := t downto 2 do { as late as possible } begin out := (dⱼ if i = m, else xᵢ₊₁,ⱼ) + yᵢ,ⱼ₊₁; xᵢⱼ := min(uᵢⱼ, out); { what capacity allows } yᵢⱼ := out − xᵢⱼ; { the rest, from j − 1 } end; xᵢ₁ := (d₁ if i = m, else xᵢ₊₁,₁) + yᵢ₂ − yᵢ₁; j := 1; while xᵢⱼ < 0 and j < t do { more on hand than needed: } begin yᵢ,ⱼ₊₁ := yᵢ,ⱼ₊₁ − xᵢⱼ; { carry it forward } xᵢ,ⱼ₊₁ := xᵢ,ⱼ₊₁ + xᵢⱼ; { and make that much less } xᵢⱼ := 0; j := j + 1; end; if xᵢⱼ < 0 or xᵢ₁ > uᵢ₁ then return infeasible; end; return x, y;end;
Model 2 implementation: Greedy as-late-as-possible procedure
# Model: greedy as-late-as-possible procedurefunctionlatestplan(d, u, y⁰, yᵗ⁺¹) m, t =size(u) x, y =zeros(m, t), zeros(m, t +1)for i in m:-1:1 on = i == m ? d : x[i+1, :] # what stage i passes on y[i, 1], y[i, t+1] = y⁰[i], yᵗ⁺¹[i] # initial, final inventoryfor j in t:-1:2# as late as possible out = on[j] + y[i, j+1] # flow out of node (i, j) x[i, j] =min(u[i, j], out) # what capacity allows y[i, j] = out - x[i, j] # the rest, from j - 1end x[i, 1] = on[1] + y[i, 2] - y[i, 1] # initial inventory first j =1while x[i, j] <0&& j < t # more on hand than needed: y[i, j+1] -= x[i, j] # carry it forward x[i, j+1] += x[i, j] # and make that much less x[i, j] =0 j +=1end (x[i, j] <0|| x[i, 1] > u[i, 1]) &&returnnothingendreturn (x = x, y = y)end
latestplan (generic function with 1 method)
Example 3: Two stages with unavailable capacity
For the product of Fig. 1, made in two stages over three months, determine the production plan that minimizes total production and inventory cost, and the cost: by hand, and with Model 1 in Julia.
A single product is produced in two stages. The planning period is three months, and the model is run each month using current demand forecasts and scheduled production capacity. In some periods, the capacity is input as zero, which indicates that it will be used for other purposes during that period. Fig. 1 shows the production-inventory planning network.
Assuming an annual inventory carrying rate of 0.3 $/$-yr, the monthly rate is h = 0.3/12. The inventory cost per ton at stage i is then the monthly rate times the cumulative cost to produce each ton from stage 1 through i:
Eq. 2 is the Units check of Lecture 1.2 at work: a carrying rate in $/$-yr becomes an inventory cost in $/ton-month only once it is divided by 12 months a year and multiplied by a cumulative cost in $/ton.
Example 3(a): By hand
Determine the plan and its cost by hand, from the capacities and flow balance.
With no inventory at stage 2 at the start or wanted at the end, and no capacity at stage 2 after the first period, all of stage 2’s production has to be in the first period. Can a ton be moved from stage 1 to stage 2 within the same period? It can: a ton can be produced at stage 1 and passed on to stage 2 in the same period, going all the way through, from raw material to finished product, within one period, and nothing forces a wait between the stages. If a stage did have to wait for the next period, the network would have a diagonal arc from stage 1 in one period to stage 2 in the next, and a process taking two months, such as some continuous processes, would have an arc spanning two periods. Networks like that very quickly get complex, but the network remains the best way to visualize them, and about the only one. The problem is also easy to make infeasible. Here all of the production has to be in the first period, where stage 1 can produce 50 tons; the total demand is less than 50 tons, but with a total demand over 50 the plan would not be feasible. That is one of the uses of these models: when demand is high, an infeasible model shows that the operation as it stands does not have the capacity, and that perhaps overtime should be scheduled. Identifying what causes the infeasibility shows what has to change to make the model feasible.
Model 2 is the flow balance of Sec. 2.2, node by node from the last, with one choice. At each node the flow out is already known: what goes on to the next stage, or to the demand at the last stage, plus the inventory carried into the next month. Production takes as much of it as the capacity allows, and the rest has to come in as inventory from the month before, which is then the inventory out of the node to its left. In month 1 the inventory coming in is the initial inventory, so production is the rest, which has to fit within the capacity. Fig. 4 draws the steps, and Table 1 lists them:
Show the code that draws this figure
# Code block 5: Model 2 by hand, a frame a node# The network with the flows found so far: its data in red, a flow still# to be found as its symbol (a production flow with its capacity), a flow# found in blue, and the node balanced at the step circled in bluefunctionplanframe!(ax, X, Y, d, u, hl; r =0.16, ls =13) m, t =size(u)hidedecorations!(ax); hidespines!(ax)P(i, j) = ((j -1) *2.1, -(i -1) *1.7) # stage i, period jarc(p, q) =arrows2d!(ax, [p[1]], [p[2]], [q[1] - p[1]], [q[2] - p[2]]; color = ink, shaftwidth =1.0, tipwidth =6, tiplength =6)tx(p, s, al) =text!(ax, p[1], p[2]; text = s, fontsize = ls, align = al)xl(i, j) =ismissing(X[i, j]) ?rich(sym("x", "$i,$j"), up(" ≤ $(u[i, j])"; color = datac)) :rich(sym("x", "$i,$j"), up(" = $(num(X[i, j]))"; color = flowc))yl(i, j) =ismissing(Y[i, j]) ? sym("y", "$i,$j") :rich(sym("y", "$i,$j"), up(" = $(num(Y[i, j]))"; color = flowc))for j in1:t # production and demand x0, y0 =P(1, j)arc((x0, y0 +0.9), (x0, y0 + r))tx((x0 +0.08, y0 +0.5), xl(1, j), (:left, :center))for i in2:m p, q =P(i -1, j), P(i, j)arc(p .- (0, r), q .+ (0, r))tx((p[1] +0.08, (p[2] + q[2]) /2), xl(i, j), (:left, :center))end xm, ym =P(m, j)arc((xm, ym - r), (xm, ym -0.7))tx((xm, ym -0.8), up("$(d[j])"; color = datac), (:center, :top))endfor i in1:m # inventory p =P(i, 1) .- (0.75, 0)arc(p, P(i, 1) .- (r, 0))tx(p .- (0.08, 0), up("$(num(Y[i, 1]))"; color = datac), (:right, :center))# each label below its arc, clear of the production label abovefor j in1:t-1 p, q =P(i, j) .+ (r, 0), P(i, j +1) .- (r, 0)arc(p, q)tx(((p[1] + q[1]) /2, p[2] -0.1), yl(i, j +1), (:center, :top))end q =P(i, t) .+ (0.75, 0)arc(P(i, t) .+ (r, 0), q)tx(q .+ (0.08, 0), up("$(num(Y[i, t+1]))"; color = datac), (:left, :center))endfor i in1:m, j in1:t # nodes, over the arcs h = (i, j) == hlpoly!(ax, Circle(Point2f(P(i, j)...), r); color =:white, strokecolor = h ? flowc : ink, strokewidth = h ? 2.6:1.4)endxlims!(ax, -1.45, (t -1) *2.1+1.4)ylims!(ax, -(m -1) *1.7-1.2, 1.0)return axend# The frames: the network as given, then one per node in the order the# procedure balances them, the last stage first and, within it, the last# period first; each fills in its node's production and inventory infunctionlatestframes(d, u, y⁰, yᵗ⁺¹) m, t =size(u) pl =latestplan(d, u, y⁰, yᵗ⁺¹) X =Matrix{Union{Missing, Float64}}(missing, m, t) Y =Matrix{Union{Missing, Float64}}(missing, m, t +1) Y[:, 1], Y[:, t+1] = y⁰, yᵗ⁺¹ nodes = [(i, j) for i inreverse(1:m) for j inreverse(1:t)] fig =Figure(size = (680, 230* (length(nodes) ÷2+1)))for k in0:length(nodes)if k ==0 hl, title = (0, 0), "(a) As given"else i, j = nodes[k] X[i, j] = pl.x[i, j] j >1&& (Y[i, j] = pl.y[i, j]) hl, title = (i, j), "($('a'+ k)) Step $k: stage $i, month $j"end ax =Axis(fig[k ÷2+1, k %2+1]; aspect =DataAspect(), title = title, titlesize =13, titlefont =:regular)planframe!(ax, X, Y, d, u, hl)endcolgap!(fig.layout, 8); rowgap!(fig.layout, 4)return figendlatestframes([20, 10, 15], [50050; 6000], [0, 0], [35, 0])
Figure 4: Ex. 3 solved by hand with Model 2, a frame per node: (a) the network as given, its data in red and each flow still to be found as its symbol, a production flow with its capacity; (b) to (g) the node balanced at each step, circled in blue, the last stage first and, within each stage, the last month first, and the production and the inventory in that it determines, in blue.
Table 1: Ex. 3 solved by hand with Model 2, a step per node: the flow out of the node, the production that takes as much of it as the capacity allows, and the inventory that comes in for the rest.
Step
Node
Flow out
Production
Inventory in
1
stage 2, month 3
d_3 + y_{2,4} = 15 + 0 = 15
x_{2,3} = \min(0, 15) = 0
y_{2,3} = 15
2
stage 2, month 2
d_2 + y_{2,3} = 10 + 15 = 25
x_{2,2} = \min(0, 25) = 0
y_{2,2} = 25
3
stage 2, month 1
d_1 + y_{2,2} = 20 + 25 = 45
x_{2,1} = 45 - 0 = 45 \le 60
y_{2,1} = 0, initial
4
stage 1, month 3
x_{2,3} + y_{1,4} = 0 + 35 = 35
x_{1,3} = \min(50, 35) = 35
y_{1,3} = 0
5
stage 1, month 2
x_{2,2} + y_{1,3} = 0 + 0 = 0
x_{1,2} = \min(0, 0) = 0
y_{1,2} = 0
6
stage 1, month 1
x_{2,1} + y_{1,2} = 45 + 0 = 45
x_{1,1} = 45 - 0 = 45 \le 50
y_{1,1} = 0, initial
By Model 2, the last stage first, stage 2 produces all three months’ demand in the first month, the only month it has capacity, x_{2,1} = 20 + 10 + 15 = 45 tons, and carries y_{2,2} = 45 - 20 = 25 tons into month 2 and y_{2,3} = 25 - 10 = 15 tons into month 3. Stage 1 feeds it in the same month, x_{1,1} = 45 tons, within its capacity of 50. Stage 1’s final inventory of 35 tons can be produced only in month 1, where 5 tons of capacity are left, or in month 3; produced in month 3, it is carried no month before the end of the horizon, so x_{1,3} = 35 tons, carried out as y_{1,4} = 35. Every other flow is zero, and the cost is
# Code block 6: Ex. 3's datau = [50050; # capacity of stage i in period j (ton)6000];d = [20, 10, 15] # demand in period j (ton)cᵖ = [200, 800] # production cost in stage i ($/ton)h =0.3/12# monthly inventory carrying rate (1/month)@show cⁱ =cumsum(cᵖ)*h # inventory cost for stage i ($/ton)y⁰ = [0, 0] # initial inventory of stage i (ton)yᵗ⁺¹ = [35, 0] # final inventory of stage i (ton)m =length(cᵖ) # number of production stages = 2t =length(d) # number of periods of production = 3LBʸ = [y⁰ zeros(m, t-1) yᵗ⁺¹] # same LB and UB to fix inventory levelUBʸ = [y⁰ fill(Inf, m, t-1) yᵗ⁺¹];
cⁱ = cumsum(cᵖ) * h = [5.0, 24.999999999999996]
The printed inventory costs are where the 5 and the 25 of Fig. 1 come from. The initial inventory is 0 at both stages and the final inventory 35 and 0, taken from the figure. The number of production stages, m, is the number of elements in the production cost, and the number of periods, t, the length of the demand vector, rather than being hard-coded, so that a change in the number of stages or periods changes them too. The lower and upper bounds on inventory, LBʸ and UBʸ, are both the initial inventory, y⁰, in their first column; between it and the last column the lower bound is zero, since inventory need not be carried, and in this problem the upper bound is infinity. Setting the lower and upper bounds to the same value clamps the inventory: it forces the variable to have that value. The final inventory, in the last column, is forced the same way, both bounds set to it. Those are all the inputs the model needs, and solving it is much like solving the minimum cost network flow problem, an application of what has already been seen.
The 25 prints as 24.999999999999996 because 0.3/12 has no exact binary representation, a rounding far too small to matter. Logjam’s snapvals, which Model 1’s implementation applies to its solution, rounds a value that close to a whole number to it: snapvals(cⁱ) is [5.0, 25.0].
Code block 7 solves Ex. 3 with the implementation of Model 1, in place of a model written out for the example, and again with Model 2’s, which needs no solver:
# Code block 7: Ex. 3 solved with prodinv, and as late as possibleTCᵒ, xᵒ, yᵒ =prodinv(cᵖ, cⁱ, d, u, LBʸ, UBʸ)@show TCᵒxˡ, yˡ =latestplan(d, u, y⁰, yᵗ⁺¹)@showsum(cᵖ .* xˡ) +sum(cⁱ .* yˡ[:, 2:end]); # its cost
1
prodinv returns a named tuple, (TC = …, x = …, y = …), and assigning it to three names unpacks its fields in order, so TCᵒ, xᵒ and yᵒ hold the total cost, the production and the inventory.
2
Model 2 finds the plan a second way, without a solver, and its cost should match TCᵒ: two routes to one number, a validation, the Triangulate check of Lecture 1.2.
Fig. 5 shows the optimal production and inventory flows (in blue) in the production-inventory planning network. Note that the flow is balanced at each node; that is, the flow into a node equals the combined production, inventory, and demand flow out of the node.
Show the code that draws this figure
# Code block 8: the optimal flows of Ex. 3 on its networkdrawpinet(cᵖ, u, cⁱ, d, y⁰, yᵗ⁺¹; X = xᵒ, Y = yᵒ)
Figure 5: The optimal plan of Ex. 3 on its network: each production and inventory flow in blue, with the costs, the initial and final inventories and the demands of Fig. 1. At every node the flow in equals the flow out.
Raw math programming output is not typically very interpretable, so part of the job of a logistics engineer, in reporting to a manager, is to add a reporting function at the end that converts the results into whatever makes sense.
The optimal production flows, with the required demands beneath them, are:
Solved with the solver’s output turned off, and checked to be optimal, the plan costs $53,175. The values are hard to read as a list, so Fig. 5 puts them in blue on the network of Fig. 1. Stage 1, with a capacity of 50 tons in the first month, produces 45 tons then and carries no inventory, and stage 2 carries 25 tons of inventory into the second month and 15 into the third. Everything balances. At stage 1 in the first month, 0 tons come in as inventory and 45 are produced, while 0 go out as inventory and 45 go on to stage 2, so the node nets to zero. At stage 2, 0 tons come in as inventory and 45 are produced, and 25 go out as inventory, so the node nets 20 tons out, the demand. Every node nets to zero except those of the last stage, which net to the demand: the plan is just a set of accounting identities that have to hold. Code blocks 9 and 10 make the values more readable, the production with the demand beneath it, and the inventory, which has four columns for three months because it runs to t + 1.
The plan is the one found by hand in Ex. 3(a), and by Model 2 in Code block 7, and so is its cost: two routes to one number, the Triangulate check of Lecture 1.2.
The same plan, at a total cost of $53,175.
Example 3(c): By hand, with capacity that binds
Ex. 3’s product, with its costs and demand, has 10 tons of inventory at stage 1 to start, none wanted at either stage at the end, and in place of those of Fig. 1 the capacities, in tons,
The steps are those of (a). This time stage 2’s capacity in month 3 is less than its flow out, so 5 of the 15 tons come in as inventory from month 2, and at stage 1 the 10 tons on hand come in to month 1, so production there is the rest. Fig. 6 draws the steps, and Table 2 lists them:
Show the code that draws this figure
# Code block 11: Ex. 3(c) by hand, a frame a nodelatestframes([20, 10, 15], [40200; 301510], [10, 0], [0, 0])
Figure 6: Ex. 3(c) solved by hand with Model 2, in the frames of Fig. 4, with the capacities of Ex. 3(c) and the 10 tons on hand at stage 1.
Table 2: Ex. 3(c) solved by hand with Model 2, a step per node.
Step
Node
Flow out
Production
Inventory in
1
stage 2, month 3
d_3 + y_{2,4} = 15 + 0 = 15
x_{2,3} = \min(10, 15) = 10
y_{2,3} = 5
2
stage 2, month 2
d_2 + y_{2,3} = 10 + 5 = 15
x_{2,2} = \min(15, 15) = 15
y_{2,2} = 0
3
stage 2, month 1
d_1 + y_{2,2} = 20 + 0 = 20
x_{2,1} = 20 - 0 = 20 \le 30
y_{2,1} = 0, initial
4
stage 1, month 3
x_{2,3} + y_{1,4} = 10 + 0 = 10
x_{1,3} = \min(0, 10) = 0
y_{1,3} = 10
5
stage 1, month 2
x_{2,2} + y_{1,3} = 15 + 10 = 25
x_{1,2} = \min(20, 25) = 20
y_{1,2} = 5
6
stage 1, month 1
x_{2,1} + y_{1,2} = 20 + 5 = 25
x_{1,1} = 25 - 10 = 15 \le 40
y_{1,1} = 10, initial
What goes in balances what goes out over the horizon, the Balance check of Lecture 1.2: the 15 + 20 + 0 tons stage 1 produces and the 10 it has on hand are the 45 tons stage 2 makes, and stage 2’s 20 + 15 + 10 = 45 tons are the 45 tons of demand. The 10 tons on hand at the start add nothing to the cost, since their cost was incurred in the previous period and the inventory sum of Eq. 1 starts at j = 2. The cost is
Produce 15 and 20 tons at stage 1 in months 1 and 2, and 20, 15 and 10 tons at stage 2 in months 1 to 3, at a total cost of $43,200.
NoteCheck Ex. 3(c) using Julia
Code block 12 solves Ex. 3(c) with Model 2’s implementation and again with Model 1’s, and the two should agree with each other and with the hand calculation, the Triangulate check of Lecture 1.2:
# Code block 12: Ex. 3(c) solved as late as possible, and with prodinvuc = [40200; # capacity of stage i in period j (ton)301510]y⁰c, yᵗ⁺¹c = [10, 0], [0, 0] # initial and final inventory (ton)xˡc, yˡc =latestplan(d, uc, y⁰c, yᵗ⁺¹c)@showsum(cᵖ .* xˡc) +sum(cⁱ .* yˡc[:, 2:end]); # its costLBʸc = [y⁰c zeros(m, t-1) yᵗ⁺¹c]UBʸc = [y⁰c fill(Inf, m, t-1) yᵗ⁺¹c]TCc, xᵒc, yᵒc =prodinv(cᵖ, cⁱ, d, uc, LBʸc, UBʸc)@show TCc;
Not every result is needed to know the plan. In Ex. 3, if the other three flows at stage 1 in month 3 are known, the only possible result for its production is the 35 tons, because every node except those of the last stage has to net to zero: whatever comes in has to go out. The only difference at the last stage is the demand that has to be satisfied.
NoteEx. 3 as a minimum cost network flow
As a minimum cost network flow, Model 3 in Lecture 4.2, the network of Fig. 1 needs a node for each arc to start from: a supply node feeding stage 1 in every period, and for each stage a node holding its initial inventory and a node receiving its final inventory. The last stage’s nodes are the demand nodes, each with the period’s demand as a negative net supply, and every other node is a transshipment node. Each production arc is bounded by its capacity, and the inventory arcs between periods are unbounded. Code block 13 repeats Lecture 4.2’s implementation of the model:
# Code block 13: list2incid and mcnf, as implemented in Lecture 4.2usingSparseArrayslist2incid(i, j, n) =sparse([i; j], [1:length(i); 1:length(i)], [ones(length(i)); -ones(length(j))], n, length(i))functionmcnf(A, c, s; u =fill(Inf, size(A, 2)), l =zeros(size(A, 2))) V, E =1:size(A, 1), 1:size(A, 2) model =Model(HiGHS.Optimizer)@variable(model, l[k] ≤ x[k ∈ E] ≤ u[k]) # (b)@objective(model, Min, sum(c[k] * x[k] for k ∈ E))@constraint(model, supply[i ∈ V; s[i] >0],sum(A[i, k] * x[k] for k ∈ E) ≤ s[i]) # (a)@constraint(model, balance[i ∈ V; s[i] ≤0],sum(A[i, k] * x[k] for k ∈ E) == s[i]) # (a)set_silent(model)optimize!(model)iftermination_status(model) == JuMP.OPTIMALreturn (TC =objective_value(model), x =snapvals(value.(x)))enderror(termination_status(model))end
mcnf (generic function with 1 method)
The initial inventory arcs cost nothing, because their cost was already incurred in the previous period; the final ones cost c_i^i, as the inventory sum of Eq. 1 runs from j = 2 to t + 1. Code block 14 builds the network’s arc list from Ex. 3’s data, production arcs first, period by period, and then each stage’s inventory arcs, solves the network, and compares its production arcs, the first mt in the list, with Ex. 3(b)’s plan:
# Code block 14: Ex. 3's time-expanded network as an arc list, solved# Node 1 supplies stage 1; node nd(i, j) is stage i in period j; then a# node per stage for its initial inventory, and one for its finalnd(i, j) =1+ (j -1)*m + iB(i) =1+ m*t + iF(i) =1+ m*t + m + in =1+ m*t +2mI, J, c, ub, lb =Int[], Int[], Float64[], Float64[], Float64[]arc!(a, b, cc, uu, ll =0.0) = (push!(I, a); push!(J, b); push!(c, cc);push!(ub, uu); push!(lb, ll))for j in1:t, i in1:m # productionarc!(i ==1 ? 1:nd(i -1, j), nd(i, j), cᵖ[i], u[i, j])endfor i in1:marc!(B(i), nd(i, 1), 0.0, y⁰[i], y⁰[i]) # initial inventoryfor j in1:t-1arc!(nd(i, j), nd(i, j +1), cⁱ[i], Inf) # inventoryendarc!(nd(i, t), F(i), cⁱ[i], yᵗ⁺¹[i], yᵗ⁺¹[i]) # final inventoryends =zeros(n) # net suppliess[1] =sum(d) +sum(yᵗ⁺¹) # enough for any planfor i in1:m s[B(i)], s[F(i)] = y⁰[i], -yᵗ⁺¹[i]endfor j in1:t s[nd(m, j)] =-d[j]end@showlength(I) # arcs@show n # nodesTCᵐ, xᵐ =mcnf(list2incid(I, J, n), c, s; u = ub, l = lb)@show TCᵐreshape(xᵐ[1:m*t], m, t) == xᵒ # the plan of Ex. 3(b)
1
Each call appends one arc to the list: start node, end node, cost, upper bound and lower bound (0 unless given). Giving the initial and final inventory arcs the same lower and upper bound fixes their flow, as LBʸ and UBʸ do in Model 1’s implementation.
length(I) = 14
n = 11
TCᵐ = 53175.0
true
The minimum cost network flow reaches the same plan at the same cost, a third route to the number found in Ex. 3(a) and (b), the Triangulate check of Lecture 1.2 again, and the demonstration that Model 1 is a special case of the minimum cost network flow problem.
The same plan, at a total cost of $53,175.
2.4 Limits of the greedy procedure
Example 4: Scheduled maintenance
A single product is produced in a two-stage production process. Stage one has a capacity to produce 40 tons of the product per month, and there is a cost of $200 for each ton of product produced. For stage two, capacity is 80 tons per month, and there is a cost of $800 per ton. Due to scheduled maintenance, no capacity will be available during month five for stage 1 and during months three and five for stage 2. Demand for the next six months is estimated to be 25, 15, 10, 50, 25, and 15 tons per month. Currently, five tons of material is in storage at stage 1 and 25 tons at stage 2; at the end of the planning period, ten tons of material should be in storage at stage 1 and 15 at stage 2. Up to 30 tons of material per month can be stored at each stage. Assuming that the inventory costs are $12 and $46 per ton-month, determine the production plan that minimizes total costs over the planning horizon and what the total production costs will be, and what the total inventory costs will be.
Solution approach: Can use the Single-Product with Variable Capacity model, Model 1. The maximum inventory storage constraint can be modeled as an upper bound on y.
Ex. 4 is how a problem would arise in practice: a description gives the demand, all of the costs and the other data, and reading it translates it into the model’s data. Here there are two stages and six periods, six months, with a demand for each. Stage 2 is shut down in month 3, and both stages in month 5, for scheduled maintenance; otherwise they run at their normal capacities of 40 and 80 tons. Instead of a carrying rate from which to calculate the inventory costs, the inventory costs c_i^i themselves are given, $12 and $46 per ton-month, along with the initial and final inventories. What is different from Ex. 3 is that the upper bound on inventory is not infinity: up to 30 tons can be stored at each stage. A limit other than infinity on how much inventory can be held between stages is not unusual: a production system does not have infinite storage, and there may be only a certain buffer for the inventory passed from one stage to the next, so a maximum inventory that can be carried is quite common. With it, the inventory’s bounds are the initial inventory, the final inventory, and an upper bound of 30 tons on every period in between, row (e) of Eq. 1.
# Code block 15: Ex. 4's datau = [40404040040; # capacity of each stage (ton)8080080080]d = [25, 15, 10, 50, 25, 15] # demand in each period (ton)# variable production cost for each stage ($/ton)cᵖ = [200, 800]cⁱ = [12, 46] # inventory cost for each stage ($/ton)y⁰ = [5, 25] # initial inventory of each stage (ton)yᵗ⁺¹ = [10, 15] # final inventory of each stage (ton)m =length(cᵖ) # number of production stagest =length(d) # number of periods of production# same LB and UB to fix initial and final inventory levelsLBʸ = [y⁰ zeros(m, t-1) yᵗ⁺¹]UBʸ = [y⁰ fill(30, m, t-1) yᵗ⁺¹]; # max inventory storage level
Negative zeros, and values that are really close to zero, are fixed by snapping each value to the closest integer when it is really close, since it is typically supposed to be an integer value, with an absolute tolerance of 10^{-8}. Logjam’s snapvals does this.
Code block 16 solves it with Model 1’s implementation, which reads every value through snapvals, as in Lecture 2.7. It then tries Model 2, whose assumption (c) the storage limit puts in doubt, and solves the problem once more without the limit:
# Code block 16: Ex. 4 solved with prodinv, and as late as possibleTCᵒ, xᵒ, yᵒ =prodinv(cᵖ, cⁱ, d, u, LBʸ, UBʸ)@show TCᵒxˡ, yˡ =latestplan(d, u, y⁰, yᵗ⁺¹)TCˡ =sum(cᵖ .* xˡ) +sum(cⁱ .* yˡ[:, 2:end]) # its cost@show TCˡ@showmaximum(yˡ[:, 2:t]) # the most it holds at a stageUBⁿ = [y⁰ fill(Inf, m, t-1) yᵗ⁺¹] # no storage limit@showprodinv(cᵖ, cⁱ, d, u, LBʸ, UBⁿ).TC;
The latest plan carries 35 tons at stage 1 into month 4, built ahead for the 75 tons stage 2 makes in month 4, that month’s demand and the next one’s, when stage 2 is down, which is more than stage 1 can make in one month; and 35 tons is more than the 30 a stage can store. Without the limit it would be optimal, at $133,900, the cost of the linear program with the limit removed; with the limit, the linear program holds part of the build-ahead at stage 2 instead, where it costs more to carry, and its plan costs $400 more. A storage limit is where the procedure stops being enough.
The optimal production flows, with the required demands beneath them, are:
The sum of the component costs exceeds the total cost from the solver because it includes the first-period (time 0) inventory costs:
TCᵖ + TCⁱ, TCᵒ # Code block 20: their sum, and the total
(135510.0, 134300.0)
Exclude first period from inventory costs:
TCⁱ =sum(cⁱ.*yᵒ[:,2:end]) # Code block 21: inventory from period 2
3300.0
The component costs now sum to the total cost from the solver:
TCᵖ + TCⁱ, TCᵒ # Code block 22: the sum, and the total
(134300.0, 134300.0)
Code blocks 19 to 22 are a sanity check, a way of making sure that a model is doing its accounting correctly, and the way the need to include the initial or the final inventory, but not both, comes to light. The total production cost and the total inventory cost are calculated from the plan, added up, and compared with the model’s objective. Summed over every inventory arc, they come to $135,510, more than the model’s $134,300, because the model excludes the first period. Recalculated from period 2, the inventory cost is $3,300 rather than $4,510, and the sum then equals the objective. For any moderately complex production-inventory problem, this bookkeeping check is worth doing almost always: whatever the model reports as its objective is recalculated from its plan, to make sure it can be duplicated. It is a great way of catching any errors that creep into the calculation, and it is the Triangulate check of Lecture 1.2, two routes to one number.
Total cost $134,300: production $131,000 and inventory $3,300.
A stage can have more than one resource. Between two nodes of the network there can then be two arcs, a multigraph, as when there are two transport options; even for a single product, one stage of production can be the transportation of the product rather than production as such, with a base transport option and a costlier overflow option. The same holds for production proper: a base machine, and, for excess production, a secondary machine that is more expensive to run but is available for overflow. That is the way an electric utility serves its demand, and it is a valuable way of looking at these problems. An electric utility typically has base load units, typically nuclear and coal units, and gas turbines, for example, that handle the peak demands. The idea is used throughout production systems as a way of handling variability in demand.
In such a multigraph the base unit typically has an upper bound, and the more expensive alternative is needed only for a demand that exceeds it. Ex. 6(d) in Lecture 4.2 is the example: the cheaper of two modes carries at most 20 units, and only what it cannot carry takes the costlier one. Overtime is such a costlier resource: it is what Ex. 3(a) suggests scheduling when a plan comes back infeasible because the operation, as it stands, does not have the capacity.
3. Capacity determined inside the model
Even where capacity is set to zero and the decision is made outside the model, it is still quite hard for a person to look at all the different trade-offs once the number of stages and periods gets long enough. Table 3 summarizes the single-product problems of this lecture by what it takes to solve each: when a simple greedy procedure solves the problem, when it is a linear program that the greedy procedure cannot solve, when it requires dynamic programming, and when it requires more. It also motivates the mixed-integer linear programming version of the problem, as a way of internalizing the need to make these decisions. In Sec. 2, the indicator that gates the capacity is provided outside the model.
Table 3: Single-product problems in this lecture: whether capacity’s availability is set outside the model or gated inside it, whether production carries a setup cost, and what solves each.
Capacity
Problem
Solved by
Set outside the model: u_{ij} given, 0 where unavailable
no setup cost, Model 2’s assumptions hold (Exs. 3, 7)
greedy procedure, Model 2; its plan is optimal, costing the same as the linear program’s
no setup cost, a storage limit binds (Ex. 4)
linear program, Model 1; the greedy plan breaks the limit
no setup cost, a base resource and a costlier one, a multigraph (Sec. 2.4)
minimum cost network flow on the multigraph, as in Lecture 4.2
Gated inside the model by indicators z_{ij}
no setup cost (Sec. 3)
linear program, Model 1: the indicators determine nothing
setup cost, one stage, no capacity limit (Sec. 3 note)
The indicator gates capacity in row (c) of Eq. 3 below, x_{ij} \le u_i z_{ij}, and the setup cost puts it in the objective, as c_i^s z_{ij}.
Two costs do not vary with the amount produced, a switchover cost and a setup cost. A switchover cost is incurred when production switches from one product to another; as long as the same product continues to be produced in the next period, there is no additional cost. A setup cost is the estimated cost of setting up to produce a product in a period. Neither is the “fixed cost” of the Location lectures, the intercept of a regression line fitted to production costs (Lecture 2.4, Sec. 8), and that fixed cost does not carry over to production-inventory planning, because production-inventory planning is primarily a tactical problem. Location is strategic: its production cost is an aggregate, and fitting a regression to capture economies of scale makes sense there. Here a period can be as short as a week or a day, too short for economies of scale, and a product made period after period would be charged the intercept again in every period. A regression’s intercept is appropriate only for strategic problems whose period is a year or more.
Most of the time, the cost of operating and producing in a period is mostly a variable cost, a constant applied to the units produced in the period. What takes the model from a linear program to a mixed-integer linear program is a setup cost in a period, assumed to be incurred in each period in which a resource produces: literally getting the machine up and running, for example, a cost independent of how much is produced but incurred every period.
In some situations, though, the setup cost is nearly everything. A truckload is one: between two points, most of its cost is independent of the load carried, up to the truck’s maximum payload, which is why Ex. 6 treats a truck as a setup.
A setup that does not have to be repeated each period is a switchover. For a single product its cost is incurred once, at the very beginning, whatever plan is made, so from a modeling point of view it can be ignored; switchover costs matter once there are several products to switch among (Lecture 4.4). Setup costs recur, each shift say, as in the validation or calibration of equipment.
A setup cost needs a variable that says whether a stage produces in a period: a binary variable z, called an indicator variable, a kind of variable used in several contexts, 1 if the stage produces in the period and 0 otherwise. Added to the objective function, it incurs the setup cost whenever the stage produces. Even with a single product, a large setup cost can change the decision: where the linear program would produce a small quantity in some month, the setup cost, independent of the quantity, can make it optimal not to produce at that stage that month, and to produce a little more earlier, carrying it as inventory, or a little later. Which periods to produce in is then determined by the model, through its indicator variables, rather than by the person running it, as the zero capacities of Sec. 2 required.
The number of truckloads in a week is not restricted to one, and typically there is no reason to restrict it: one dedicated truck making round trips might take all week, but in practice a third-party trucking company can send one, two or a dozen trucks at any time. So the number of trucks is an integer variable rather than a binary one, and the setup cost is incurred for however many truckloads are sent.
Model 3 generalizes Model 1 with the number of setups at each stage in each period. A machine is set up at most once in a period, so for a machine the number is the indicator; for trucks it is a count.
minimize: total production, inventoryand setup cost over the planning horizon
solve for:
(a) production at each stage in each period;
(b) inventory at each stage carried into each period; (c) number of setups at each stage in each period.
subject to:
(a) flow balance at each stage that feeds the next: what is produced or carried in is either passed on to the next stage or carried out;
(b) flow balance at the last stage, which feeds demand: what is produced or carried in either meets the period’s demand or is carried out;
(c) capacity: production at each stage in each period is no more than its capacityper setup times its number of setups;
(d) initial and final inventory: inventory at each stage at the start and at the end of the horizon is given;
(e) storage: inventory at each stage between periods is no more than its storage capacity.
assumptions:
(a) production cost and inventory cost per ton at each stage are constant;
(b) demand in each period is given;
(c) production can pass through every stage within the period it starts in; (d) setup cost is incurred for each setup at a stage, whatever the setup then produces; (e) inventory carries each ton’s share of its setups, taken at a full setup’s capacity.
The other symbols are those of Eq. 1. Both lines labeled (d) in Eq. 3 are row (d), as in Eq. 1.
What to use for the inventory cost per period is not obvious once production has a setup cost as well as a cost proportional to the tons produced. The setup cost cannot simply be left out of it, because it can be the most significant cost relative to inventory: in some situations a setup can cost more than the value of the product, being essentially the main production activity, and dropping it would make such situations grossly inaccurate. The setup cost is allocated to the inventory cost at a stage by assumption (e) of Model 3. Each ton is charged a share of every setup it has passed through, the setup cost divided by the capacity per setup, added to the production cost of each stage up to the stage it is at before the carrying rate is applied:
Using the capacity u_i in the denominator is conservative with respect to the cost allocation, since most stages probably operate below capacity: a setup is then shared by fewer tons than u_i, and Eq. 4 charges each ton the smallest share it could carry. It slightly understates the inventory cost, but it includes some part of the setup cost incurred every time a stage produces in a period.
Some inputs are “lumpy,” that is, the cost of utilizing the input is the same over a range of outputs (e.g., the labor costs to set up a machine are the same whether one or 100 parts are produced).
Constraint (c) represents available production capacity and uses the setups z to indicate whether production occurs. The capacity is assumed to be available in discrete units that are either at 100% of the unit’s capacity or not available at all. Integer decision variables are used to make these decisions. This internal capacity decision requires the model to be solved as a MILP.
Producing anything therefore requires a setup, and a setup is paid for in the objective. For a machine, \bar{z}_{ij} is 1 in a period in which it is available and 0 in one in which it is down; for trucks, it is as many as could ever be needed. With no setup costs and machines only, Model 3 is Model 1 with capacities u_{ij} = u_i\bar{z}_{ij}.
Time spent setting up a machine when the machine is stopped and can not be operating is termed the machine’s internal setup time; in contrast, the machine’s external setup time is the time required to prepare to produce a particular part on the machine that does not require the machine to be stopped (e.g., fixture preparation). A machine’s “down time” represents the percentage of time the machine is not available for production due to repair or planned maintenance; it does not include internal setup time because the machine is still available for production during this time.
Setup costs are also where Model 2 stops applying,2 and they make the problem far harder to solve.
A mixed-integer linear program like this one would have been hard to solve in the 1980s for a production problem of any realistic size; it was not until the early to mid 1990s that such problems became feasible to solve.
Even a problem like this one with a setup cost, solved in a minute, can change: in a mixed-integer program, changing one number from 100 to 300 can leave the problem feasible but take it from solving in one minute to three hours with the solver, since it is highly nonlinear.
Instead of capacities changed by hand period by period, each stage has a single capacity for all periods, and the indicator variables act like a switch, or a valve, turning that capacity off or on; otherwise the model is essentially the same as before, and with the setup cost allocated to the inventory cost by Eq. 4, everything else carries through as before.
The implementation is almost the same as before. The setups are added to the objective, multiplied by the setup cost, and the capacity, which could be a bound in the definition of the production variable while it was given to the model, becomes a constraint of its own, since it is now determined in the model: the capacity per setup times the setups turns the capacity off or on. In the solution the setups should match the production, zero where nothing is produced and positive where something is: the production gives the amount, and the setups say whether the stage produces in the period. The accounting check is done as before. As a linear program, it was up to the person running the model to work out in which periods to produce; with the setup cost in the model, the model makes that decision itself, and that is the gateway to much more complex decision-making.
Model 3 implementation: Production-inventory with setups
# Model: production-inventory planning with setup costsfunctionprodinv_setup(cᵖ, cˢ, cⁱ, d, u, z̄, LBʸ, UBʸ) m, t =size(z̄) model =Model(HiGHS.Optimizer)@variable(model, x[i=1:m, j=1:t] ≥0)@variable(model, 0≤ z[i=1:m, j=1:t] ≤ z̄[i,j], Int)@variable(model, LBʸ[i,j] ≤ y[i=1:m, j=1:t+1] ≤ UBʸ[i,j])@objective(model, Min,sum(cᵖ[i]*x[i,j] + cˢ[i]*z[i,j] for i=1:m, j=1:t) +sum(cⁱ[i]*y[i,j] for i=1:m, j=2:t+1))@constraint(model, flowbal[i=1:m-1, j=1:t], # (a)-x[i,j] + x[i+1,j] - y[i,j] + y[i,j+1] ==0)@constraint(model, laststage[j=1:t], # (b)-x[m,j] - y[m,j] + y[m,j+1] ==-d[j])@constraint(model, capacity[i=1:m, j=1:t], # (c) x[i,j] ≤ u[i]*z[i,j])set_silent(model)optimize!(model)iftermination_status(model) == JuMP.OPTIMALreturn (TC =objective_value(model), x =snapvals(value.(x)), y =snapvals(value.(y)), z =Int.(snapvals(value.(z))))enderror(termination_status(model))end
1
The stages and periods are read from the most setups, z̄[i, j], a matrix laid out like the capacities u[i, j] of Model 1’s implementation; the capacity per setup u is now a vector, one value for each stage.
prodinv_setup (generic function with 1 method)
Without the setup cost. Take c_i^s = 0 in Eq. 3. The indicator still gates the capacity, row (c), but nothing in the objective makes the model turn it off, so wherever a stage may be set up its capacity is there at no cost: the model is Model 1 with u_{ij} = u_i\bar{z}_{ij}, and which periods are available is determined outside it, as in Sec. 2.
For this problem, a single product, there is no easy way of making the decision inside the model without something like a setup cost. Additional constraints could be added to the problem to make the gate determine something, but a model is kept as simple as possible, with elements added only when they are absolutely needed, and adding to the model just to make the decision is extra complexity rather than the problem’s natural state. The setup cost is the natural element: it is what makes producing in a period cost something, and so what makes in which periods to produce a decision for the model.
With several products, the subject of Lecture 4.4, a gate alone does determine something. When a resource is shared exclusively among products, its capacity is turned on or off by the same production indicator, and a linking constraint lets only one product’s indicator be 1 at a stage in a period: one product or the other, not both. The objective function does not include any additional terms. The linking constraint is an additional constraint, with nothing added to the objective function, and everything else flows through as before.
Example 5: Setup costs in manufacturing
A single product is produced in a two-stage production process. Stage one has the capacity to produce 40 tons of the product per month, and there is a setup cost of $6,000 per month and a variable cost of $200 for each ton of product produced. For stage two, capacity is 80 tons per month, and there is a setup cost of $18,000 and a variable cost of $800. The setup costs are only incurred when any amount of production occurs at that stage for that month. Demand for the next six months is estimated to be 25, 15, 10, 50, 25, and 15 tons per month. Currently, five tons of material are in storage at stage one and 25 tons at stage two; at the end of the planning period, ten tons of material should be in storage at stage one and 15 at stage two. Up to 30 tons of material per month can be stored at each stage. Assuming that the annual inventory carrying rate is 0.4 $/$-yr, determine the production plan that minimizes total costs over the planning horizon and what those total costs will be.
Solution approach:Model 3 extends Model 1 with a production indicator binary variable to turn capacity on or off at each stage for each period. One modeling issue is determining inventory cost per ton if the production cost includes a setup cost. An approximation is to allocate the setup cost based on capacity, Eq. 4, each stage’s setup cost c_i^s divided by its capacity u_i.
The maximum inventory storage constraint can be modeled as an upper bound on y.
# Code block 23: Ex. 5's datau = [40, 80] # capacity of each stage (ton)d = [25, 15, 10, 50, 25, 15] # demand in each period (ton)# variable production cost for each stage ($/ton)cᵖ = [200, 800]cˢ = [6000, 18000] # setup cost for each stage ($/month)h =0.4/12# monthly inventory carrying rate (1/month)# inventory cost for each stage ($/ton)@show cⁱ =cumsum(cᵖ + cˢ./u)*hy⁰ = [5, 25] # initial inventory of each stage (ton)yᵗ⁺¹ = [10, 15] # final inventory of each stage (ton)m =length(u) # number of production stagest =length(d) # number of periods of production# same LB and UB to fix initial and final inventory levelsLBʸ = [y⁰ zeros(m, t-1) yᵗ⁺¹]UBʸ = [y⁰ fill(30, m, t-1) yᵗ⁺¹]; # max inventory storage level
The optimal setups, the production indicators, are:
zᵒ # Code block 27: setups
2×6 Matrix{Int64}:
0 1 1 1 0 1
0 1 0 1 0 1
Calculate total variable production, inventory, and setup cost components:
# Code block 28: the cost components, inventory from period 2TCᵖ, TCⁱ, TCˢ =sum(cᵖ.*xᵒ), sum(cⁱ.*yᵒ[:,2:end]), sum(cˢ.*zᵒ)
(131000.0, 3275.0000000000005, 78000)
The component costs sum to the total cost from the solver:
TCᵖ + TCⁱ + TCˢ, TCᵒ # Code block 29: their sum, and the total
(212275.0, 212275.0)
Incurring the setup cost removes decisions to produce just a little, say 5% of capacity in a month: in some months nothing is produced, to save the setup cost, and inventory is carried instead, which is cheaper than producing every month.
Here stage 1 runs in four months and stage 2 in three, each skipping month 1, where the initial inventory covers the demand, and month 5, where production in month 4 is carried instead.
Total cost $212,275: production $131,000, inventory $3,275, and setups $78,000.
Example 5(b): Cost of a given plan by hand
Using some different data as noted below (otherwise, use the data given in Ex. 5) and only a three-month planning horizon, determine the total variable production, inventory, and setup cost over the horizon. The forecasted demand is 22, 3, and 33 tons per month. The product loses five percent of its value after one month. Each month, 20, 4, 40 and 25, 0, 37 tons, respectively, are produced at each stage. Currently, five tons of material is in storage at stage one and none at stage two; at the end of the planning period, seven tons of material should be in storage at stage one and four at stage two.
The production values, the inputs to the network, are given, so the missing data can be inferred using the flow balance equations, as in Exs. 1 and 2, and the costs then computed from them.
Going month by month, the inventory carried out of each node is what came in less what went on:
The final inventories come out at the 7 and 4 tons required, the Balance check of Lecture 1.2. The product losing five percent of its value in a month is a carrying rate of h = 0.05 per month, and with the setup share of assumption (e) the inventory costs are c_1^i = 0.05\,(200 + 6{,}000/40) = 17.50 and c_2^i = 0.05\,(200 + 6{,}000/40 + 800 + 18{,}000/80) = 68.75 $/ton-month. Stage 1 produces in all three months and stage 2 in months 1 and 3, so
Code block 30 computes the same three costs from the plan, and checks them against the answers published with the problem, the Landmark check of Lecture 1.2:
# Code block 30: Ex. 5(b)'s plan, costedX = [20440; 25037] # production at stage i in month jdb = [22, 3, 33] # demand in month jY =zeros(2, 4) # inventory, by flow balanceY[:, 1] = [5, 0]for j in1:3 Y[1, j+1] = Y[1, j] + X[1, j] - X[2, j] Y[2, j+1] = Y[2, j] + X[2, j] - db[j]end@assert Y[:, 4] == [7, 4] # the final inventories requiredcⁱb =cumsum(cᵖ + cˢ./u)*0.05# monthly carrying rate 0.05Z = X .>0# a setup wherever there is productionVPC, IC, SC =sum(cᵖ.*X), sum(cⁱb.*Y[:, 2:end]), sum(cˢ.*Z)@assertall((VPC, IC, SC) .≈ (62400, 673.75, 54000)) # published@show VPC@show IC@show SC;
VPC = 62400
IC = 673.75
SC = 54000
Variable production cost $62,400, inventory cost $673.75, and setup cost $54,000, for a total of $117,073.75.
What drives the decision to produce in a period is whether its setup cost is worth incurring, and what counts as a setup cost varies. In baking, firing up the oven costs about the same for one loaf of bread as for a hundred, so that cost can be treated as a setup cost, while the variable cost is mostly the ingredients: the more bread, the more flour. Whether to fire up the oven on a given day or sell bread from the day before is the same trade-off as in Ex. 5, although bread goes stale, so the oven is probably fired up every day; other resources would go unused in some periods. Transportation is the same phenomenon, set up differently. Thought of in more general terms, transportation is a kind of production, the production of a movement, and a truckload is at the other extreme from the oven. Moving a tractor-trailer from one point to another costs very nearly the same with one box in it as completely full, since most of the cost is the truck itself and its driver, who is paid either way; only the fuel varies slightly, mostly with wind resistance. Easily 90 or 95 percent of the cost is independent of the quantity in the truck,3 so as a good approximation it is all setup cost.
The Transport lectures give a transport stage its costs. The truckload charge can be viewed as the setup cost, with effectively almost no variable cost, since the charge is essentially for the truck and trailer. Ex. 6 is the transport counterpart of Ex. 5’s manufacturing.
Example 6: Truckload as a setup cost
The product of Ex. 1 in Lecture 3.3, 20 tons a year of 40-lb, 9-ft³ cartons valued at $25,000 a ton and carried at 0.3 $/$-yr, is shipped by truckload 532 miles from Raleigh (27606) to Gainesville (32606) at January 2018 rates. Plan its shipments week by week over a year as Model 3, with a demand of 20/52 tons a week at Gainesville and no inventory at either end of the year: (a) determine the plan; and (b) determine it again at a value of $85,000 a ton.
The plan has two stages. Stage 1 is the plant, which can produce any amount in a week at the product’s value, $25,000 a ton, with no setup cost. Stage 2 is the shipment by truckload, set up once for each truck sent, at one truck’s charge, with no cost per ton and a capacity per setup of the maximum payload; any number of trucks can be sent in a week. The inventory at Gainesville carries the product’s value plus the truckload charge per ton, assumption (e) of Model 3, which is Lecture 3.3’s allocated full-truckload rate, Ex. 1(h).
The cost of having a truck is found from the truck rate, and that cost becomes the setup cost of the production-inventory problem: every time a truck comes, its cost is incurred, the same kind of cost as a setup.
Code block 31 repeats the data of Lecture 3.3, Ex. 1, and the charge for one truck:
# Code block 31: Lecture 3.3's Ex. 1, and one truck's chargeuwt, ucu =40, 9# lb and ft^3 per cartons = uwt/ucu # density (lb/ft^3)Kwt, Kcu =25, 2750# trailer limits (ton, ft^3)@show qmax =maxpayld(s; Kwt = Kwt, Kcu = Kcu) # tonf =20# ton/yrd =532# mi, Raleigh to GainesvilleppiTL, ppiLTL =131.0, 177.4# PPIs for Jan 2018rTL =2.00ppiTL/102.7# $/micTL =charge_tl(qmax, s, d; r = rTL, Kwt = Kwt, Kcu = Kcu, ppi = ppiTL) # $/truck@show cTL;
The truckload’s setup cost is one truck’s charge, $1,357.2, with no cost per ton, and its capacity per setup is the maximum payload, 6.11 tons.
Example 6(a): Weekly plan
Determine the weekly plan for the year.
Code block 32 builds the year’s data for Model 3 from the truck’s charge and payload, as a function of the product’s value v, so that part (b) can change it, and solves it:
# Code block 32: the year planned week by week, at product value vusingStatisticsfunctionweekly(v) weeks, h =52, 0.3# weeks in the year; carrying rate dw =fill(f/weeks, weeks) # demand at Gainesville (ton/week) cᵖ = [v, 0.0] # plant; truck ($/ton) cˢ = [0.0, cTL] # setup costs ($/setup) u = [f, qmax] # capacity per setup (ton) z̄ = [ones(Int, 1, weeks); # most setups in a week: the plant,fill(10, 1, weeks)] # and trucks, more than ever needed cⁱ = h/weeks .* [v, v + cTL/qmax] # inventory cost ($/ton-week) LBʸ =zeros(2, weeks +1) # no inventory at either end UBʸ = [zeros(2) fill(Inf, 2, weeks -1) zeros(2)]returnprodinv_setup(cᵖ, cˢ, cⁱ, dw, u, z̄, LBʸ, UBʸ)endr25 =weekly(25000)tl25 =findall(>(0), r25.z[2, :]) # weeks a truck is sent@showlength(tl25)@showmean(r25.x[2, tl25]) # tons per truckmean(diff(tl25)) # weeks between trucks
The plan sends 7 trucks in the year, about 2.86 tons each, one every 7.5 weeks on average. Lecture 3.3 found the same shipment by an entirely different route, minTLC, its total logistics cost minimized over a continuous shipment size.
As Lecture 3.3 put it, most companies more realistically do batch production: a customer’s product is produced over a very short time, so there is essentially no inventory at the origin. With consumption still constant at the destination, \alpha = 0 + 1/2 = 1/2, where before, with inventory at both ends, it was 1.
Code block 33 asks minTLC for the shipment with \alpha = 1/2:
# Code block 33: Lecture 3.3's minTLC on the same shipmenttr = (r = rTL, Kwt = Kwt, Kcu = Kcu, ppi = ppiTL)res25 =minTLC((f = f, s = s, a =1/2, v =25000, h =0.3, d = d), tr, ppiLTL)@show res25.qᵒ # ton@show52res25.qᵒ/f # weeks between shipmentsres25.isLTL
res25.qᵒ = 2.690423228838891
(52 * res25.qᵒ) / f = 6.995100394981117
false
minTLC ships 2.69 tons by truck every 7.0 weeks. Two models, a mixed-integer linear program in whole weeks from Networks and a continuous minimization from Transport, reach nearly the same plan, the Triangulate check of Lecture 1.2; the weekly plan’s trucks are a little larger because it can ship only whole weeks’ demand, and must end the year with nothing on hand.
# Code block 34: the plan and minTLC at a value of $85,000 a tonr85 =weekly(85000)tl85 =findall(>(0), r85.z[2, :])@showlength(tl85)@showmean(r85.x[2, tl85])@showmean(diff(tl85))res85 =minTLC((f = f, s = s, a =1/2, v =85000, h =0.3, d = d), tr, ppiLTL)@show res85.qᵒ52res85.qᵒ/f
At the higher value, inventory costs more to carry, so the shipments should be smaller and more frequent, and they are: 13 trucks of about 1.54 tons, one every 4.0 weeks, against minTLC’s 1.46 tons every 3.8 weeks.
Getting the general trends right is always important with computer work, because it is easy to put numbers in and get crazy results; a mental model of the trends to expect is what flags a result that has gone wrong.
An answer that moves the right way when an input is raised is the Nudge check of Lecture 1.2.
13 trucks of about 1.54 tons each.
4. Forecasting demand
Every model so far has taken the demand in each period, d_j, as given. In practice, the initial periods of a planning horizon are probably occupied primarily by orders already booked: the orders already received are used to populate demand. Toward the end of the horizon the orders thin out, and rather than demand decreasing over the horizon, new orders will be received. That is where traditional forecasting is applied, based on some type of recent history, to forecast the likely orders, especially for products that have seasonality or are not expected to be constant.
The demand over the planning horizon typically comes from a corporate information technology environment, a master production schedule or a material requirements planning (MRP) or enterprise resource planning (ERP) system, which holds the best information available. That information is a mix: for some periods the demand is just a forecast, and for others there are firm orders, customers promised delivery six months out, say. In some businesses the demand going forward is the aggregate of the orders actually received, and little forecasting is done at all: everything to be produced over the next year is produced against firm orders in hand. In a company, the demand comes down a planning hierarchy. Sales and operations planning sets a plan for each family of products, typically monthly and a year or more out; master scheduling turns it into a master production schedule (MPS) for each product, which says what will be made in each period of the planning horizon; and MRP uses the bills of material to work out every component the MPS needs and when to order it. Going down the hierarchy, the detail increases and the horizon shrinks.4 An enterprise resource planning (ERP) system is a software system that controls the entire operations of a firm. The item, carrier, and customer master files are maintained by an ERP system and are used as a common data source for orders.
Of the three scope levels of logistics decisions in Lecture 1.1, Sec. 2, strategic (years), tactical (weeks to a year) and operational (minutes to a week), the production-inventory models of Secs. 2 and 3 are tactical. The demand is a mix of firm and forecasted orders: further out, more of it is forecast, and the more immediate periods have more actual booked orders. Planning guides draw the same shape as a curve, the cumulative orders on hand displacing the forecast the closer a period comes to the present.5 An order that has been booked was part of the forecast, so it cannot simply be added to it. Planning systems let orders consume the forecast: within a short window starting at the present, the demand time fence, only the orders booked count; beyond it, one of several common rules uses the larger of the forecast and the orders already booked.6Fig. 7 draws both: the orders on hand displacing the forecast as a period comes closer, and the rule on either side of the fence.
Show the code that draws this figure
# Code block 35: orders on hand consuming the forecastlet F = [0.85, 1.10, 0.80, 0.95, 1.00, 0.70, 0.95, 1.20, 1.05, 0.80] n, fence =length(F), 2# periods ahead, demand time fence ordc, fcc = colorant"#1f5fa8", colorant"#b9cde5" grey = colorant"#bdbdbd" s = [max(0.0, 1- (ℓ -1)/6)^1.4 for ℓ in1:n] # share usually booked b = s .* F # orders on hand b[4] =1.18# booked beyond its forecast fig =Figure(size = (680, 350)) ax =Axis(fig[1, 1]; xlabel ="Periods ahead", ylabel ="Demand", xticks =1:n)hideydecorations!(ax; label =false)for ℓ in1:npoly!(ax, Rect2f(ℓ -0.4, 0, 0.8, b[ℓ]); color = ordc)if b[ℓ] < F[ℓ] && ℓ <= fence # forecast set asidepoly!(ax, Rect2f(ℓ -0.4, b[ℓ], 0.8, F[ℓ] - b[ℓ]); color =:white, strokecolor = grey, strokewidth =1)elseif b[ℓ] < F[ℓ] # forecast not yet bookedpoly!(ax, Rect2f(ℓ -0.4, b[ℓ], 0.8, F[ℓ] - b[ℓ]); color = fcc)elseif b[ℓ] > F[ℓ] # the forecast, below the orderslines!(ax, [ℓ -0.4, ℓ +0.4], [F[ℓ], F[ℓ]]; color =:white, linestyle =:dash, linewidth =1.6)endendtext!(ax, 4, F[4] +0.03; text ="forecast", fontsize =11, color =:white, align = (:center, :bottom))text!(ax, 4, b[4] +0.04; fontsize =11, align = (:center, :bottom), text ="orders exceed\nthe forecast:\nplan on the orders")vlines!(ax, [fence +0.5]; color = ink, linestyle =:dash, linewidth =1.2)text!(ax, fence +0.6, 1.85; text ="demand time fence", fontsize =12, align = (:left, :bottom))text!(ax, (fence +1)/2, 1.38; text ="orders only", fontsize =12, align = (:center, :bottom))text!(ax, 7.6, 1.38; text ="the larger of\nforecast and orders", fontsize =12, align = (:center, :bottom))arrows2d!(ax, [9.6], [0.42], [-5.8], [0]; color = ink, shaftwidth =1.2, tipwidth =8, tiplength =8)text!(ax, 7.0, 0.46; fontsize =11, align = (:center, :bottom), text ="as a period comes closer, orders replace its forecast") elem = [PolyElement(color = ordc), PolyElement(color = fcc),PolyElement(color =:white, strokecolor = grey, strokewidth =1)]Legend(fig[2, 1], elem, ["Orders on hand", "Forecast", "Forecast set aside"]; orientation =:horizontal, framevisible =false, labelsize =12)xlims!(ax, 0.4, n +0.6); ylims!(ax, 0, 2.05); hidespines!(ax, :t, :r) figend
Figure 7: Schematic of the orders on hand consuming the forecast, for a demand that varies from period to period. The nearer a period, the more of its demand is already booked (dark blue) and the less is forecast (light blue). Within the demand time fence only the orders count, and the forecast not yet booked is set aside; beyond it the plan uses the larger of the forecast and the orders, so in period 4, where the orders already exceed the forecast (dashed), it plans on the orders.
The more immediate periods typically have better forecasts, and the further into the future, the less accurate the forecasts become. A forecast is always wrong, so it should come with a measure of its error, measured at the lead at which the forecast is acted on, since the error grows with the lead; and simple methods are hard to beat.78 How long a horizon should be, the subject of Sec. 5, depends on that growth of the error with the lead.
5. Planning horizon
How is the planning horizon determined? It has to be long enough to incorporate any decision that affects the production of the product over time. If something plays itself out in two or three weeks, a horizon of a month is more than adequate, whereas a production process that has to be shut down for month-long maintenance may need a horizon of two or three years, to incorporate several of its maintenance downtimes into the plan.
Not every product needs a long horizon. One that takes three or four days to make, from materials in hand or easy to get, on resources available most of the month, needs no plan of several months: a week ahead is enough, and that is production scheduling, not a problem for these models. The models are for decisions that take months to play out. One is maintenance: resources that each have to be down for two months every three years need to be taken down at the best times, never in the peak season and staggered so that together they disrupt the least, given the best estimate of demand. Another is lead time: with a material that takes four or five months to arrive once ordered, production waits until the material is in hand rather than holding work in process for months.
Oil refineries, for example, plan month by month but typically about 36 months ahead, a window within which each refinery is taken down for one, two or three months of scheduled maintenance.9 Given an estimate of demand, the models work out when to take each resource down so that several are not down at once, which makes them a planning tool for a complex organization.
The length of the horizon should include all significant planning decisions, and one of the main ones is scheduled maintenance. In one study of power plants in Saudi Arabia, maintenance was done in the winter, when demand was lower than in the summer, and the plan looked about three years ahead to determine which plants would be down for maintenance, and when, while still meeting demand.
A typical number of periods is anywhere from six to 36; with monthly periods, three years of monthly data is common. With weekly periods, a plan might look 26 weeks ahead.
The research literature calls a horizon long enough that nothing beyond it can change the first period’s decision a forecast horizon.1011 Four questions follow, in turn: how a plan uses its horizon, how long the horizon must be, why it should not be longer, and how a choice of horizon is checked. Fig. 8 draws the first two for a resource with an outage ahead, its materials ordered two periods before the period that uses them.
Show the code that draws this figure
# Code block 36: rolling plans ahead of an outage, with a lead timelet T, δ =12, 2# periods, lead time cap, o =0.7, 9:10# capacity, periods of the outage d = [0.50, 0.33, 0.46, 0.37, 0.47, 0.41,0.42, 0.41, 0.41, 0.42, 0.43, 0.44] # demand planc, ordc = colorant"#b9cde5", colorant"#1f5fa8" idlec, invc = colorant"#9ccf9c", colorant"#e08a1e" grey = colorant"#bdbdbd" nb, sp =0, 0.0# count back from the outagewhile sp <sum(d[o]) nb +=1 sp += cap - d[o[1] - nb]end b0 = o[1] - nb # first period of the build-ahead a0 = b0 - δ # its materials are ordered H = δ + nb +length(o) rows = T - H +1# plans whose window fits fig =Figure(size = (680, 640))# (a) a plan made in each period: its first δ periods already# ordered, the period δ ahead ordered now, the rest planned ax1 =Axis(fig[1, 1]; xticks =1:T, ylabel ="Plan made in period", yticks =1:rows, yreversed =true, titlesize =13, titlefont =:regular, title =rich("(a) Rolling: plan ",it("H"), up(" = $H"), " periods ahead; ","order for the period a lead time away"))hidespines!(ax1, :t, :r)vspan!(ax1, [b0 -0.5], [o[1] -0.5]; color = (idlec, 0.35))vspan!(ax1, [o[1] -0.5], [o[end] +0.5]; color = (datac, 0.18))for k in1:rows J = k:min(k + H -1, T)for j in Jif j < k + δ # already orderedpoly!(ax1, Rect2f(j -0.43, k -0.30, 0.86, 0.60); color =:white, strokecolor = ordc, strokewidth =1.2)elsepoly!(ax1, Rect2f(j -0.45, k -0.32, 0.9, 0.64); color = j == k + δ ? ordc : planc, strokecolor =:white)endend# inventory carried in at the window's start, out at its endscatter!(ax1, [k -0.5], [k]; marker =:rtriangle, markersize =9, color = invc)last(J) == k + H -1&&scatter!(ax1, [last(J) +0.5], [k]; marker =:rtriangle, markersize =9, color = invc)endpoly!(ax1, Rect2f(a0 -0.5, a0 -0.42, H, 0.84); color =:transparent, strokecolor = ink, strokewidth =1.4)text!(ax1, 0.6, a0 +1.6; fontsize =11, align = (:left, :center), text ="the outlined plan, made\nin period $a0, "*"orders for\nperiod $b0, the first of the\n"*"build-ahead, so it has\nto see past the outage")xlims!(ax1, 0.4, T +0.6); ylims!(ax1, rows +0.6, 0.4) elem = [PolyElement(color =:white, strokecolor = ordc, strokewidth =1.2),PolyElement(color = ordc), PolyElement(color = planc),MarkerElement(marker =:rtriangle, color = invc, markersize =9)]Legend(fig[2, 1], elem, ["Already ordered", "Ordered now", "Planned, not yet ordered","Inventory carried in at the start and out at the end"]; orientation =:horizontal, nbanks =2, framevisible =false, labelsize =12)# (b) the horizon of the outlined plan: from its order to the end# of the outage, each period's materials ordered δ periods ahead ax2 =Axis(fig[3, 1]; xlabel ="Period", xticks =1:T, titlesize =13, titlefont =:regular, title ="(b) How long: from the first build-ahead "*"order to the end of the outage")hideydecorations!(ax2); hidespines!(ax2, :t, :r, :l)for j in1:Tif b0 <= j < o[1] # idle capacity builds inventorypoly!(ax2, Rect2f(j -0.45, d[j], 0.9, cap - d[j]); color = idlec)endpoly!(ax2, Rect2f(j -0.45, 0, 0.9, d[j]); color = j in o ? datac : grey)endlines!(ax2, [0.5, o[1] -0.5, o[1] -0.5, o[end] +0.5, o[end] +0.5, T +0.5], [cap, cap, 0, 0, cap, cap]; color = ink, linewidth =1.5)for k in a0:(o[1] -1- δ) # one order a period, δ ahead xs =range(k, k + δ; length =40)lines!(ax2, xs, cap +0.04.+0.13.*sin.(π.* (xs .- k) ./ δ); color = ordc, linewidth =1.4)scatter!(ax2, [k], [cap +0.04]; markersize =5, color = ordc)scatter!(ax2, [k + δ], [cap +0.04]; marker =:dtriangle, markersize =8, color = ordc)endtext!(ax2, a0 + δ/2, cap +0.19; text =it("δ"), fontsize =13, color = ordc, align = (:center, :bottom))text!(ax2, a0 -0.6, cap +0.21; fontsize =11, color = ordc, align = (:right, :center), text ="one order each\nperiod, a lead\ntime ahead")text!(ax2, T +0.4, cap +0.02; text ="capacity", fontsize =12, align = (:right, :bottom))text!(ax2, 0.6, 0.05; text ="demand", fontsize =12, align = (:left, :bottom))text!(ax2, (b0 + o[1] -1)/2, cap -0.03; fontsize =11, align = (:center, :top), text ="idle capacity\nbuilds inventory")text!(ax2, (o[1] + o[end])/2, d[o[end]] +0.03; text ="outage", fontsize =12, color = datac, align = (:center, :bottom))bracket!(ax2, a0 -0.45, cap +0.33, o[end] +0.45, cap +0.33; text =rich(it("H"), up(" ≥ "), it("δ"), up(" + "),"build-ahead", up(" + "), "outage"), fontsize =12, style =:square, color = ink, textcolor = ink, textoffset =18)xlims!(ax2, 0.4, T +0.6); ylims!(ax2, 0, 1.4)linkxaxes!(ax1, ax2)rowsize!(fig.layout, 1, Relative(0.5)) figend
Figure 8: Schematic of rolling plans ahead of an outage, with materials ordered two periods before the period that uses them. (a) A plan is made in each period over a horizon of H = 7 periods: its first two periods were ordered by the plans before it (outlined), it orders the materials for the period two ahead (dark blue), and the rest is planned but not yet ordered (light blue); each plan starts from the inventory carried in and ends with the inventory carried out (orange). (b) The idle capacity of periods 6 to 8 (green) builds the inventory that covers the demand of the outage in periods 9 and 10 (red). Each period orders for the period two ahead, so the plan made in period 4, outlined in (a), is the earliest that has to see to the end of the outage, and H \ge \delta + build-ahead + outage.
A plan is acted on only in its first period. For Ex. 3’s plan, the manager is told to produce 45 tons at each stage in period 1. The later periods are not acted on, but they are part of the planning: how much to produce now depends on how much will be carried over as inventory and how much produced later. A month later the model is run again, and the second period becomes the first. Even if no demand has changed, a planned maintenance may have moved and freed capacity, which changes the plan from the inventory carried over; and each new run adds an estimate of the demand one more period out, which can change what is done now. Problems of realistic size give almost no insight just by looking at them.
The work flow of a planning model in use, as Fig. 9 draws it, runs in five steps:
Reports are run on a regular period-to-period, rolling-horizon basis as part of normal operations management. The model is only changed when the logistics network changes.
Show the code that draws this figure
# Code block 37: the logistics software stacklet boxc, edgec = colorant"#dff2da", colorant"#3a8a4f"# name, what it is, centre (x, y); every box 6.65 wide and 2.5 high B = [("Data", "csv, Excel, etc.", (4.0, -1.3)), ("Report", "GUI, web, etc.", (11.3, -1.3)), ("Scripting", "Python, Matlab, Julia", (7.62, -5.92)), ("User Library", "in script language", (16.15, -5.92)), ("MIP Solver", "Gurobi, Cplex, etc.", (4.0, -10.55)), ("Standard Library", "in compiled C, Java", (11.3, -10.55))] w, h, rr =6.65, 2.5, 0.9# width, height, cornerfunctionrbox(c) # a rounded box outline x0, y0 = c[1] - w/2+ rr, c[2] - h/2+ rr x1, y1 = c[1] + w/2- rr, c[2] + h/2- rrarc(xc, yc, a) = [Point2f(xc +rr*cos(θ), yc +rr*sin(θ)) for θ inrange(a, a +π/2; length =12)] cs = [(x1, y0, -π/2), (x1, y1, 0), (x0, y1, π/2), (x0, y0, π)]returnreduce(vcat, [arc(k...) for k in cs]) # corner by cornerend fig =Figure(size = (680, 430)) ax =Axis(fig[1, 1]; aspect =DataAspect())hidedecorations!(ax); hidespines!(ax)# what commercial software packages, dashedlines!(ax, [0.05, 19.95, 19.95, 0.05, 0.05], [-3.65, -3.65, -12.35, -12.35, -3.65]; color = ink, linestyle =:dash, linewidth =1.5)text!(ax, 17.45, -9.55; text ="Commercial\nSoftware\n(LLamasoft, etc.)", align = (:center, :center), justification =:center, fontsize =15, color = ink)# each arrow from box edge to box edge; the last three run both waysfor (p, q, both) in [((4.0, -2.55), (6.15, -4.67), false), # data in ((9.25, -4.67), (11.25, -2.55), false), # report ((10.945, -5.92), (12.825, -5.92), true), ((6.1, -7.17), (4.1, -9.3), true), ((9.3, -7.17), (11.25, -9.3), true)]for (a, b) in (both ? [(p, q), (q, p)] : [(p, q)])arrows2d!(ax, [a[1]], [a[2]], [b[1] - a[1]], [b[2] - a[2]]; color = datac, shaftwidth =2.2, tipwidth =11, tiplength =11)endendfor (name, what, c) in B # boxes, over the arrowspoly!(ax, rbox(c); color = boxc, strokecolor = edgec, strokewidth =2)text!(ax, c[1], c[2] +0.45; text = name, fontsize =17, font =:bold, color = ink, align = (:center, :center))text!(ax, c[1], c[2] -0.55; text ="($what)", fontsize =14, color = ink, align = (:center, :center))endxlims!(ax, -0.2, 20.2); ylims!(ax, -12.6, 0.0) figend
Figure 9: The logistics software stack, the work flow of a planning model in use: the data go into a script, which calls a MIP solver and its libraries and turns the output into a report. Commercial software, such as LLamasoft’s, packages everything inside the dashed line.
With the models of this lecture done, the overall work flow can be seen better. The first step is getting the data, and for a model run every week or every month, getting them is worth automating: fed in from whatever system the company uses, its forecast of demand, say, and whatever other information is needed, so that they are taken into the model on a regular basis without much work. A scripting language, Julia, Matlab or Python, is a good way of getting the data: it takes them in whatever form they come, from a database or elsewhere, and converts them into what the mathematical programming model needs.
The script then calls the solver, and the last piece is reporting the results in a user-friendly way. With 20 plants, not every plant manager needs the entire output of the run: each is sent what that plant should produce over the next several periods, the output focused on each entity in the logistics network. Data, model, solver, output and report is the basic flow.
A plan’s final inventory stands in for every period beyond its horizon, and with nothing required at its end, a plan too short to see a coming outage will not build for it. The cumulative lead time \delta is the longest chain of purchasing and production lead times.
Each period’s materials are ordered \delta periods ahead of the period that uses them, so the first \delta periods of every plan are already committed by orders placed in the periods before, and the decision a plan made now acts on is the order for period \delta + 1; with the materials on hand, \delta = 0.12 A plan therefore looks at least \delta periods ahead, and further when the order placed now depends on what comes after the period it supplies. For an outage, a count estimates how much further. Scheduled maintenance is one of the important reasons to plan several periods ahead: the plan has to make sure that enough inventory is on hand to cover the maintenance.
An outage with no other source. This count is the machine-hours accounting of Lecture 1.3, Sec. 5, done one period at a time. A resource at utilization u (Eq. 5 in Lecture 1.3) has idle capacity of 1 - u of what it can produce in a period. Lecture 1.3 warned that u < 1 holds only on average, and a scheduled outage is a burst of demand over capacity that is known in advance and too long to average away. The plan builds the inventory for it on the idle capacity of the periods before it: counting back from the outage, each period’s idle capacity is added until the total covers the outage’s demand, and the horizon has to span those periods and the outage. The plan that has to see the outage is the one that orders the materials for the first period of the build-ahead, \delta periods before it, the plan outlined in Fig. 8; each later period’s materials are ordered by a later plan, which sees the outage as well. At a steady utilization u the count can be done in closed form. Measure demand and inventory in periods of the resource’s production, so that each period’s demand is u and its idle capacity 1 - u. An outage of \tau periods then has a demand of \tau u, all of it to be built ahead, and b periods of build-ahead supply b\,(1 - u) of it, so the number of periods of build-ahead satisfies:
b\,(1 - u) = \tau u, \qquad b = \frac{\tau u}{1 - u}
The plan that has to see the outage spans, in turn, the \delta periods already committed at its front, the b periods of build-ahead and the \tau periods of the outage, the three spans Fig. 8(b) brackets, so its horizon is at least their sum. With b substituted, and since \tau u + \tau (1 - u) = \tau, the sum is Eq. 5:
=utilization of the resource in the periods before the outage.
By Eq. 5, with materials on hand (\delta = 0), a month-long outage at 90% utilization needs a horizon of 10 months, and at 95%, 20; a lead time adds to each. The horizon grows as 1/(1-u), the factor by which cycle time grows in Lecture 1.3’s VUT equation (Eq. 6 in Lecture 1.3): a busy resource has little idle capacity to build ahead with, which is one reason a month of maintenance can call for years of planning; the other is that the plan has to span the maintenance cycle, so that it can determine when each plant is down.
The lead time is counted once in Eq. 5, not once for each period of the build-ahead, because the materials for a build-ahead are not all ordered at one time and left to sit for two or three periods; they are ordered each period, so a two-period lead time overlaps each period of the build-ahead.
The plan is re-solved every period, so the numbers get less certain the further out they are, and are more subject to change. Each run can also change what was planned for the coming periods, which unsettles everyone who works to the plan; companies damp this by freezing the nearest periods of the plan, at some cost.13
A longer horizon can also break the model: one that solves well can go to running for 80 hours without stopping when one more variable is added, or when the horizon goes from 12 periods to 18. These models assume the demand will be satisfied exactly as planned; safety stock applies only to the current period.
Solve with H and again with a longer one: if neither the first period’s decision nor the cost of the plan acted on changes, H was long enough for these data, the Nudge check of Lecture 1.2. As Sec. 1 said, more than one horizon can be defended; what makes the choice right is its reason.
Example 7: Horizon of the month-1 plan
Ex. 3’s product, with Ex. 3’s costs, has a forecast of demand and a schedule of capacities for the next six months, the three of Ex. 3 and three more, starting and ending with no inventory. The capacities of its two stages and the demand in the six months are, in tons,
The six months are the data available, not the horizon. Determine by hand how many of them a plan made in month 1 has to look ahead, and the plan over all six months.
This is a two-stage, six-month example. It is the same as Ex. 3, except that the figure of Ex. 3’s network showed only the first three periods; otherwise the data are the same.
Starting and ending with no inventory, as this example and Ex. 6 do, is rare in practice, since these are typically ongoing operations; the zeros are a simplification.
Example 7(a): By hand
Determine the horizon from the capacities and the demand, before any plan is solved.
Stage 2 has no capacity in months 2 and 3, an outage of two months, so their demand, 10 and 15 tons, has to be made in month 1. In the terms of Eq. 5’s count, with the materials on hand (\delta = 0), stage 2’s idle capacity in month 1, 60 − 20 = 40 tons, covers the outage’s 25 tons, so the build-ahead is one month, and the plan made in month 1 has to span it and the outage: H \ge 0 + 1 + 2 = 3 months. That is an estimate committed before any plan is solved, the Prior check of Lecture 1.2.
Example 7(b): In Julia
Determine the plan over all six months.
Code block 38 enters the six months, Ex. 3’s three and the next three, and solves them with Model 1:
# Code block 38: the six-month planu6 = [500500050; # capacity of stage i in period j (ton)6000606060]d6 = [20, 10, 15, 20, 15, 25] # demand in period j (ton)cᵖ = [200, 800] # production cost in stage i ($/ton)cⁱ =cumsum(cᵖ)*0.3/12# inventory cost for stage i ($/ton)LBʸ =zeros(2, 7) # no inventory at the start or endUBʸ = [zeros(2) fill(Inf, 2, 5) zeros(2)]TC6, x6, y6 =prodinv(cᵖ, cⁱ, d6, u6, LBʸ, UBʸ)@show TC6prt(y6; rows = ["1", "2"], cols =string.(1:7), row_title ="Stage") # inventory into each month
The plan carries 35 tons of stage-1 inventory out of month 3, because stage 1 has no capacity in months 4 and 5 and stage 2 needs it to meet their demand. That is where Ex. 3’s final inventory of 35 tons at stage 1 came from: Ex. 3 is the first three months of this plan, and its final inventory stands in for the months beyond its horizon. The plan also agrees with the count of Ex. 7(a): stage 2 makes 45 tons in month 1, the demand of months 1 to 3, and carries 25 and 15 tons of it into months 2 and 3.
NoteCheck the horizon using Julia
Code block 39 checks the count of Ex. 7(a). It solves the plan made in month 1 with Model 1 over each horizon H from one month to six, with nothing required at its end, and keeps each plan’s month-1 production, the decision that is acted on:
# Code block 39: the plan made in month 1, over horizons of 1 to 6 monthsx₁ =zeros(2, 6) # month-1 production, horizon Hfor H in1:6 LBʸ =zeros(2, H +1) # no inventory at the start or end UBʸ = [zeros(2) fill(Inf, 2, H -1) zeros(2)] x₁[:, H] =prodinv(cᵖ, cⁱ, d6[1:H], u6[:, 1:H], LBʸ, UBʸ).x[:, 1]endprt(x₁; rows = ["1", "2"], cols = ["H = $H" for H in1:6], row_title ="Stage") # month 1's production, by H
Stage H = 1 H = 2 H = 3 H = 4 H = 5 H = 6
───────────────────────────────────────────────
1 20 30 45 45 45 45
2 20 30 45 45 45 45
The month-1 decision changes until H = 3 and not after. With one month, stage 2 makes only month 1’s 20 tons, and with two, 30, so the plan made in month 2, with no capacity at stage 2 in months 2 and 3, cannot satisfy their demand; from three months on, every plan makes 45 tons at stage 2 in month 1, whatever it sees beyond. So 3 months is long enough for these data, as the count found, and a longer horizon leaves the first month’s decision unchanged: the check of a horizon stated above, the Nudge check of Lecture 1.2.
A plan made in month 1 has to look 3 months ahead; the plan over all six months costs $106,250.
Endnotes
Its objective and its assumption (c), that production started in a period is completed in it, are those of the single-stage model in A. Kefeli, R. Uzsoy, Y. Fathi and M. Kay, “Using a mathematical programming model to examine the marginal price of capacitated resources,” International Journal of Production Economics 131(1), 2011, 383-391, which points, for production that takes longer than a period, to S. T. Hackman and R. C. Leachman, “A general framework for modeling production,” Management Science 35(4), 1989, 478-495.↩︎
For one stage with a setup cost and no limit on its capacity, the dynamic program of Wagner and Whitin finds the plan without a solver, since each setup makes exactly the demand of the periods up to the next setup: H. M. Wagner and T. M. Whitin, “Dynamic version of the economic lot size model,” Management Science 5(1), 1958, 89-96. With capacities, as in Ex. 5, the problem is hard in general: M. Florian, J. K. Lenstra and A. H. G. Rinnooy Kan, “Deterministic production planning: algorithms and complexity,” Management Science 26(7), 1980, 669-679.↩︎
In 2024 fuel was about 21 percent of the cost of operating a truck, $0.48 of $2.26 a mile, and only part of the fuel varies with the weight carried, so well over 80 percent of a truckload’s cost does not vary with the quantity in it: American Transportation Research Institute, An Analysis of the Operational Costs of Trucking: 2025 Update, July 2025, truckingresearch.org, accessed 2026-10-16.↩︎
APICS, CPIM Learning System, Module 1, excerpts published by ASCM: 2020 edition, pp. 1-58 to 1-66, quoting the APICS Dictionary, 16th ed., for the planning hierarchy (learningsystem.ascm.org/wp-content/uploads/2019/11/CPIM_LS_2020_Excerpt.pdf); 2019 edition, pp. 1-207 to 1-211, for orders displacing forecasts and the time fences (learningsystem.ascm.org/wp-content/uploads/2018/11/APICS_CPIM_2019_Excerpt.pdf); accessed 2026-10-16.↩︎
APICS, CPIM Learning System, Module 1, excerpts published by ASCM: 2020 edition, pp. 1-58 to 1-66, quoting the APICS Dictionary, 16th ed., for the planning hierarchy (learningsystem.ascm.org/wp-content/uploads/2019/11/CPIM_LS_2020_Excerpt.pdf); 2019 edition, pp. 1-207 to 1-211, for orders displacing forecasts and the time fences (learningsystem.ascm.org/wp-content/uploads/2018/11/APICS_CPIM_2019_Excerpt.pdf); accessed 2026-10-16.↩︎
APICS, CPIM Learning System, Module 1, excerpts published by ASCM: 2020 edition, pp. 1-58 to 1-66, quoting the APICS Dictionary, 16th ed., for the planning hierarchy (learningsystem.ascm.org/wp-content/uploads/2019/11/CPIM_LS_2020_Excerpt.pdf); 2019 edition, pp. 1-207 to 1-211, for orders displacing forecasts and the time fences (learningsystem.ascm.org/wp-content/uploads/2018/11/APICS_CPIM_2019_Excerpt.pdf); accessed 2026-10-16.↩︎
S. Chopra and P. Meindl, Supply Chain Management: Strategy, Planning, and Operation, Ch. 7, “Demand Forecasting in a Supply Chain,” the characteristics of forecasts.↩︎
R. J. Hyndman and G. Athanasopoulos, Forecasting: Principles and Practice, 3rd ed., OTexts, 2021, Secs. 5.5 and 8.3, otexts.com/fpp3, accessed 2026-10-16. Free online, and the place to read further on forecasting.↩︎
Most refineries are on a three- to five-year cycle of these maintenance shutdowns, called turnarounds, each planned years ahead and scheduled where possible in spring or fall, when the demand for gasoline is lower: American Fuel & Petrochemical Manufacturers, “Refinery turnarounds 101: What are turnarounds and why do we need them?”, October 17, 2023, afpm.org, accessed 2026-10-16.↩︎
A. Federgruen and M. Tzur, “Minimal forecast horizons and a new planning procedure for the general dynamic lot sizing model: Nervousness revisited,” Operations Research 42(3), 1994, 456-468.↩︎
C.-Y. Lee and E. V. Denardo, “Rolling planning horizons: error bounds for the dynamic lot size model,” Mathematics of Operations Research 11(3), 1986, 423, doi:10.1287/moor.11.3.423, bound the cost of planning over a finite horizon by the value of perfect information about the data beyond it.↩︎
APICS, CPIM Learning System, Module 1, excerpts published by ASCM: 2020 edition, pp. 1-58 to 1-66, quoting the APICS Dictionary, 16th ed., for the planning hierarchy (learningsystem.ascm.org/wp-content/uploads/2019/11/CPIM_LS_2020_Excerpt.pdf); 2019 edition, pp. 1-207 to 1-211, for orders displacing forecasts and the time fences (learningsystem.ascm.org/wp-content/uploads/2018/11/APICS_CPIM_2019_Excerpt.pdf); accessed 2026-10-16.↩︎
A. Federgruen and M. Tzur, “Minimal forecast horizons and a new planning procedure for the general dynamic lot sizing model: Nervousness revisited,” Operations Research 42(3), 1994, 456-468.↩︎
