MentDB WeakMQL Recipes
← innov-ai.com All profiles Portal
Profiles / Sciences
OR

Operations researcher

Optimisation with constraints, Markov chains for maintenance, newsvendor inventory: solve and plot.

10 recipes on this page · open in the interactive view

Recipes
  1. Product mix: linear program and its graph
  2. Availability of a machine with a Markov chain
  3. Newsvendor: optimal order under normal demand
  4. Call centre staffing with Erlang C
  5. PERT: probability to meet the contract deadline
  6. EOQ, safety stock and reorder point
  7. Trend and Q4 peak forecast by QR least squares
  8. Weibull bearing: age replacement policy
  9. Leontief input-output model of three sectors
  10. Discrete newsvendor from an empirical demand
executed on MentDB

Product mix: linear program and its graph

Maximise the profit 40x + 30y of two products under machine, labour and market constraints, find the optimal vertex as the intersection of the two saturated constraints, and draw the feasible region with the iso-profit line.

Pipeline
01 · optimiseMaximise with constraintsnum maximize→02 · vertexIntersection of the active constraintsnum solve_system→03 · constraintsConstraint lines and iso-profitplot function→04 · marketVertical line x = 40plot xy→05 · optimumMark the optimal vertexplot points
Script
admin
#1. Product mix: profit 40x + 30y, machine 2x + y <= 100, labour x + y <= 80, market x <= 40;
-> "[opt]" (num maximize "40*x+30*y" "{\"x\": 1, \"y\": 1}" "{\"constraints\": [\"2*x + y <= 100\", \"x + y <= 80\", \"x <= 40\"], \"bounds\": {\"x\": [0, null], \"y\": [0, null]}}");
-> "[z]" (json load select "rx_lp" [opt] /value);
#2. The optimum is the vertex where machine and labour are both saturated;
-> "[v]" (num solve_system "[\"2*x + y = 100\", \"x + y = 80\"]" "{\"x\": 10, \"y\": 10}");
json load "rx_vx" [v];
#3. Draw the constraints, the market limit, the iso-profit line and the optimal vertex;
plot function (concat "[\"100-2*x\", \"80-x\", \"(" [z] "-40*x)/30\"]") "x" 0 40;
plot add (plot xy "[40, 40]" "[0, 80]");
plot add (plot points (concat "[" (json select "rx_vx" /solution/x) "]") (concat "[" (json select "rx_vx" /solution/y) "]"))
Result
mentdb
The chart is shown in the Plot tab of the editor.
Go further
  • Read the shadow prices: increase the machine capacity by 1 and run num maximize again
  • Solve the vertex exactly with cas solve
executed on MentDB

Availability of a machine with a Markov chain

A machine is up, degraded or down (then repaired) from one day to the next. Compute the stationary law, check it with P^50, then the availability, the mean daily margin and the mean number of days before a breakdown (fundamental matrix).

Pipeline
01 · systemP transpose - I with a row of onesmatrix sub→02 · stationarySolve for the stationary lawmatrix solve→03 · checkRow of P^50matrix pow→04 · absorptionN = (I - Q)^-1matrix inverse→05 · marginExpected margin per dayvector dot
Script
admin
#1. Daily transitions of a machine: up, degraded, down - repaired;
-> "[P]" "[[0.90, 0.08, 0.02], [0.00, 0.85, 0.15], [0.60, 0.00, 0.40]]";
#2. Stationary law: pi P = pi and sum pi = 1, system Pt - I with the last row replaced by ones;
-> "[M]" (matrix sub (matrix transpose [P]) (matrix identity 3));
-> "[M]" (matrix set (matrix set (matrix set [M] 3 1 1) 3 2 1) 3 3 1);
-> "[pi]" (matrix solve [M] "[0, 0, 1]");
#3. Check: any row of P^50 converges to pi;
-> "[P50]" (matrix row (matrix pow [P] 50) 1);
#4. Mean days before a breakdown: fundamental matrix N = I - Q inverse on the up and degraded states;
-> "[Q]" "[[0.90, 0.08], [0.00, 0.85]]";
-> "[N]" (matrix inverse (matrix sub (matrix identity 2) [Q]));
-> "[T]" (matrix mul [N] "[1, 1]");
concat "stationary " [pi] ", P^50 row 1 " [P50] ", availability " (math round (math calc (concat "1-" (matrix get (concat "[" [pi] "]") 1 3))) 4) ", margin per day " (math round (vector dot [pi] "[1000, 600, -400]") 2) " EUR, days to breakdown from up / degraded " [T]
Result
mentdb
stationary [0.5882352941176472,0.31372549019607837,0.0980392156862745], P^50 row 1 [0.5882354296305502,0.313725390547185,0.09803917982226582], availability 0.902, margin per day 737.25 EUR, days to breakdown from up / degraded [15.333333333333336,6.666666666666666]
Go further
  • Compare two maintenance policies by changing the second row of P
  • Look at the speed of convergence with the second eigenvalue of matrix eigen
executed on MentDB

Newsvendor: optimal order under normal demand

A seasonal product sells 12, costs 5 and is salvaged at 1, with a normal demand (500, 120). Find the order quantity from the critical ratio with erf, confirm it by minimising the expected overage and underage cost, and compute the expected leftover and profit.

Pipeline
01 · critical ratioCR = cu / (cu + co)math calc→02 · quantilePhi(z) = CR with erfnum solve→03 · costMinimise the expected costnum minimize→04 · checkDerivative zero at q*num derivative→05 · leftoverIntegral of (q - x) f(x)num integrate
Script
admin
#1. Newsvendor: demand normal mean 500 sd 120, price 12, cost 5, salvage 1;
-> "[mu]" 500; -> "[s]" 120; -> "[p]" 12; -> "[c]" 5; -> "[v]" 1;
#2. Critical ratio CR = cu / cu + co, and q* with Phi of z = CR, Phi written with erf;
-> "[cr]" (math calc "([p]-[c])/([p]-[v])");
-> "[q]" (num solve (concat "0.5*(1+erf((q-" [mu] ")/(" [s] "*sqrt(2))))-" [cr]) "q" 0 1000);
#3. Expected cost of overage and underage with the normal loss function;
-> "[z]" (concat "((q-" [mu] ")/" [s] ")");
-> "[L]" (concat "(exp(-" [z] "^2/2)/sqrt(2*pi)-" [z] "*0.5*erfc(" [z] "/sqrt(2)))");
-> "[cost]" (concat "(" [c] "-" [v] ")*(q-" [mu] "+" [s] "*" [L] ")+(" [p] "-" [c] ")*" [s] "*" [L]);
-> "[best]" (num minimize [cost] "q" 0 1000);
#4. Checks: derivative zero at q*, expected leftover by integration;
-> "[slope]" (num derivative [cost] "q" [q]);
-> "[left]" (num integrate (concat "(" [q] "-x)*exp(-((x-" [mu] ")/" [s] ")^2/2)/(" [s] "*sqrt(2*pi))") "x" "-inf" [q]);
concat "CR = " (math round [cr] 4) ", q* = " (math round [q] 1) ", minimize: " [best] ", slope at q* = " (math round [slope] 8) ", expected leftover " (math round [left] 1) " units, expected profit " (math round (math calc (concat "([p]-[c])*" [mu] "-" (json load select "rx_nv" [best] /value))) 2)
Result
mentdb
CR = 0.6364, q* = 541.9, minimize: {"x": 541.8506834620454, "value": 495.5326352140132}, slope at q* = 0, expected leftover 71.7 units, expected profit 3004.47
Go further
executed on MentDB

Call centre staffing with Erlang C

A call centre receives 120 calls per hour with an average handling time of 4 minutes. The script increases the number of agents until the 80/20 service level is met, with Erlang B taken from the Poisson law and Erlang C, then gives the average speed of answer and the occupancy.

Pipeline
01 · trafficOffered load A = lambda.AHTmath calc→02 · erlang BB = Poisson pmf over cdf at cstat dist_pdf→03 · erlang CP(wait) and service level, loop on cstat dist_cdf→04 · kpiASA and occupancymath calc
Script
admin
#1. Call centre: 120 calls per hour, average handling time 240 s, offered traffic A in Erlangs;
-> "[lam]" 120; -> "[aht]" 240; -> "[A]" (math calc "[lam]*[aht]/3600");
#2. Target 80/20: 80 % of the calls answered within 20 s. Smallest number of agents c;
-> "[c]" (math ceiling [A]); -> "[sl]" 0;
while (< [sl] 0.8) {
	++ "[c]";
	#Erlang B from the Poisson law, then Erlang C;
	-> "[B]" (/ (stat dist_pdf "poisson" (concat "[" [A] "]") [c]) (stat dist_cdf "poisson" (concat "[" [A] "]") [c]));
	-> "[C]" (math calc "[c]*[B]/([c]-[A]*(1-[B]))");
	-> "[sl]" (math calc "1-[C]*exp(-([c]-[A])*20/[aht])");
};
#3. Average speed of answer and occupancy;
-> "[asa]" (math calc "[C]*[aht]/([c]-[A])");
concat "A = " [A] " Erlangs, agents = " [c] ", P(wait) = " (math round [C] 4) ", service level = " (math round (* 100 [sl]) 1) " % in 20 s, ASA = " (math round [asa] 1) " s, occupancy = " (math round (math calc "100*[A]/[c]") 1) " %"
Result
mentdb
A = 8 Erlangs, agents = 11, P(wait) = 0.245, service level = 80.9 % in 20 s, ASA = 19.6 s, occupancy = 72.7 %
Go further
  • Add a shrinkage of 30 % for breaks and training: scheduled agents = c / 0.7 with math ceiling
  • Check the M/M/c result with a discrete-event simulation drawing times with stat dist_random
executed on MentDB

PERT: probability to meet the contract deadline

The four tasks of a critical path have optimistic, most likely and pessimistic durations. Beta-PERT gives each expected duration as a 1-4-1 weighted mean and its variance; the normal approximation gives the probability to finish within 40 days and the 95 % date.

Pipeline
01 · PERT mean(a + 4m + b)/6stat weighted_mean→02 · varianceSum of ((b - a)/6)^2stat sum_sq→03 · probabilityP(T <= 40 d)stat dist_cdf→04 · quantileDate with 95 % confidencestat dist_inv
Script
admin
#1. PERT: optimistic, most likely, pessimistic durations in days of the critical path tasks;
-> "[w]" "[1, 4, 1]";
-> "[d1]" (stat weighted_mean "[4, 6, 10]" [w]);
-> "[d2]" (stat weighted_mean "[8, 12, 22]" [w]);
-> "[d3]" (stat weighted_mean "[3, 5, 9]" [w]);
-> "[d4]" (stat weighted_mean "[6, 9, 18]" [w]);
#2. Beta-PERT variances sigma^2 = b-a over 6 squared, summed along the path;
-> "[var]" (stat sum_sq (concat "[" (math calc "(10-4)/6") ", " (math calc "(22-8)/6") ", " (math calc "(9-3)/6") ", " (math calc "(18-6)/6") "]"));
-> "[T]" (math calc "[d1]+[d2]+[d3]+[d4]");
#3. Probability to finish within the 40-day contract, and the date with 95 % confidence;
-> "[law]" (concat "[" [T] ", " (math sqrt [var]) "]");
-> "[p40]" (stat dist_cdf "normal" [law] 40);
-> "[d95]" (stat dist_inv "normal" [law] 0.95);
concat "durations " (math round [d1] 2) ", " (math round [d2] 2) ", " (math round [d3] 2) ", " (math round [d4] 2) " d; critical path = " (math round [T] 2) " d, sigma = " (math round (math sqrt [var]) 2) " d; P(T <= 40 d) = " (math round (* 100 [p40]) 1) " %; 95 % date = " (math round [d95] 1) " d"
Result
mentdb
durations 6.33, 13, 5.33, 10 d; critical path = 34.67 d, sigma = 3.38 d; P(T <= 40 d) = 94.3 %; 95 % date = 40.2 d
Go further
  • Check the near-critical paths too: a parallel path with a larger variance can be the real risk
  • Simulate triangular durations with stat dist_random and compare the histogram with plot histogram
executed on MentDB

EOQ, safety stock and reorder point

A spare part sells 12000 units a year. The script computes the economic order quantity, then the safety stock and reorder point for a 97.5 % cycle service level, getting z from the inverse error function and checking it with erf.

Pipeline
01 · eoqQ* = sqrt(2DS/H)math calc→02 · zz = sqrt 2 erf_inv(2p - 1)math erf_inv→03 · checkPhi(z) with erfmath erf→04 · roundOrder quantitymath rint
Script
admin
#1. Spare part: 12000 units a year, 150 EUR per order, holding cost 2.40 EUR per unit and year;
-> "[D]" 12000; -> "[S]" 150; -> "[H]" 2.4;
#2. Economic order quantity Q* = sqrt 2DS/H and its yearly cost;
-> "[Q]" (math calc "sqrt(2*[D]*[S]/[H])");
-> "[TC]" (math calc "[D]/[Q]*[S]+[Q]/2*[H]");
#3. Service level 97.5 % on a 2-week lead time, weekly demand sigma 40: z from the inverse error function;
-> "[z]" (math calc (concat "sqrt(2)*" (math erf_inv (math calc "2*0.975-1"))));
-> "[SS]" (math calc "[z]*40*sqrt(2)");
-> "[ROP]" (math calc "[D]/52*2+[SS]");
#4. Check: probability of no stock-out during the lead time, Phi z = 1/2 1 + erf z/sqrt 2;
-> "[p]" (math calc (concat "0.5*(1+" (math erf (math calc "[z]/sqrt(2)")) ")"));
concat "Q* = " (math rint [Q]) " units (" (math round (math calc "[D]/[Q]") 1) " orders a year, cost " (math round [TC] 2) " EUR); z = " (math round [z] 4) ", safety stock " (math ceil [SS]) ", reorder point " (math ceil [ROP]) " units; cycle service level check " (math round [p] 4)
Result
mentdb
Q* = 1225 units (9.8 orders a year, cost 2939.39 EUR); z = 1.96, safety stock 111, reorder point 573 units; cycle service level check 0.975
Go further
  • Add quantity discounts and compare the total cost of each price break
  • Use the fill rate with the loss function instead of the cycle service level
executed on MentDB

Trend and Q4 peak forecast by QR least squares

Eight quarters of demand show a trend and a Q4 peak. The script fits y = a + b t + c Q4 with a QR decomposition (no normal equations), reads the residual sum of squares from Q^T y and forecasts the next Q4.

Pipeline
01 · decomposeX = Q Rmatrix qr→02 · projectQ^T ymatrix transpose→03 · solveR1 beta = (Q^T y)1matrix solve→04 · forecastNext Q4vector dot
Script
admin
#1. Quarterly demand of a product over 2 years, model y = a + b t + c Q4 with a Q4 peak dummy;
-> "[y]" "[410, 432, 445, 560, 470, 488, 502, 618]";
-> "[X]" "[[1, 1, 0], [1, 2, 0], [1, 3, 0], [1, 4, 1], [1, 5, 0], [1, 6, 0], [1, 7, 0], [1, 8, 1]]";
#2. QR decomposition of the design matrix, stable least squares without the normal equations;
json load "rx_p4b_qr" (matrix qr [X]);
-> "[Qty]" (matrix mul (matrix transpose (json select "rx_p4b_qr" /Q)) [y]);
-> "[R]" (concat "[" (matrix row (json select "rx_p4b_qr" /R) 1) "," (matrix row (json select "rx_p4b_qr" /R) 2) "," (matrix row (json select "rx_p4b_qr" /R) 3) "]");
#3. Back substitution R1 beta = first 3 terms of Q^T y, the other terms give the residual sum of squares;
-> "[beta]" (matrix solve [R] (concat "[" (atom get (string replace (string replace [Qty] "[" "") "]" "") 1 ",") "," (atom get (string replace (string replace [Qty] "[" "") "]" "") 2 ",") "," (atom get (string replace (string replace [Qty] "[" "") "]" "") 3 ",") "]"));
-> "[rss]" (math calc (concat (vector dot [Qty] [Qty]) "-" (vector dot (matrix mul [R] [beta]) (matrix mul [R] [beta]))));
json unload "rx_p4b_qr";
#4. Forecast of the next Q4, t = 12;
concat "a, b, c = " [beta] "; RSS = " (math round [rss] 2) "; forecast Q4 of year 3: " (math round (vector dot [beta] "[1, 12, 1]") 0) " units"
Result
mentdb
a, b, c = [399.0555555555555,14.694444444444438,101.77777777777776]; RSS = 37.47; forecast Q4 of year 3: 677 units
Go further
  • Add one dummy per quarter and compare the RSS
  • Check the coefficients with matrix solve on X directly, which also uses QR
executed on MentDB

Weibull bearing: age replacement policy

A pump bearing follows a Weibull law (beta 2.5, eta 1200 h). The script gives the MTTF with the Gamma function and scans the preventive replacement age that minimises the cost per hour against run-to-failure.

Pipeline
01 · mttfeta Gamma(1 + 1/beta)math gamma→02 · cycleExpected cycle lengthnum integrate→03 · costCost rate C(T)math calc→04 · compareSaving against run-to-failuremath round
Script
admin
#1. Pump bearing, Weibull life beta = 2.5, eta = 1200 h, preventive change 400 EUR, failure 3000 EUR;
-> "[b]" 2.5; -> "[eta]" 1200;
#2. Mean time to failure eta Gamma 1+1/beta;
-> "[mttf]" (math calc (concat [eta] "*" (math gamma (math calc "1+1/[b]"))));
#3. Age replacement: cost per hour C T = Cp R T + Cf F T over the integral of R from 0 to T;
-> "[best]" "";
-> "[cmin]" 1000000;
for (-> "[T]" 200) (<= [T] 1400) (-> "[T]" (+ [T] 50)) {
	-> "[den]" (num integrate "exp(-(t/1200)^2.5)" "t" 0 [T]);
	-> "[cT]" (math calc (concat "(400*exp(-(" [T] "/1200)^2.5)+3000*(1-exp(-(" [T] "/1200)^2.5)))/" [den]));
	if (< [cT] [cmin]) { -> "[cmin]" [cT]; -> "[best]" [T]; };
};
#4. Compare with run-to-failure, cost 3000 / MTTF;
concat "MTTF = " (math round [mttf] 1) " h; best preventive interval " [best] " h at " (math round [cmin] 4) " EUR/h, against " (math round (math calc "3000/[mttf]") 4) " EUR/h if run to failure: saving " (math round (math calc "(1-[cmin]*[mttf]/3000)*100") 1) " %"
Result
mentdb
MTTF = 1064.7 h; best preventive interval 500 h at 1.3945 EUR/h, against 2.8177 EUR/h if run to failure: saving 50.5 %
Go further
  • Refine the optimum with num minimize around the best grid point
  • Add the downtime cost and an inspection interval
executed on MentDB

Leontief input-output model of three sectors

With the technical coefficients of agriculture, industry and services, the script builds the Leontief matrix I - A, reads its determinant from the LU factors, solves the total output for the final demand and gives the output multipliers.

Pipeline
01 · leontiefI - Amatrix sub→02 · factorLU decompositionmatrix lu→03 · output(I - A) x = dmatrix solve→04 · multipliersColumn sums of the inversematrix inverse
Script
admin
#1. Input-output table of 3 sectors: agriculture, industry, services, technical coefficients A;
-> "[A]" "[[0.15, 0.10, 0.05], [0.25, 0.35, 0.15], [0.10, 0.20, 0.25]]";
-> "[d]" "[120, 300, 450]";
#2. Leontief matrix I - A and its LU factors;
-> "[M]" (matrix sub (matrix identity 3) [A]);
json load "rx_p4b_lu" (matrix lu [M]);
-> "[U]" (json select "rx_p4b_lu" /U);
-> "[det]" (math calc (concat (matrix get [U] 1 1) "*" (matrix get [U] 2 2) "*" (matrix get [U] 3 3)));
json unload "rx_p4b_lu";
#3. Total output needed for the final demand: I - A x = d;
-> "[x]" (matrix solve [M] [d]);
#4. Output multipliers: column sums of the Leontief inverse;
-> "[Li]" (matrix inverse [M]);
-> "[mult]" (matrix row (matrix mul "[[1, 1, 1]]" [Li]) 1);
concat "det(I-A) = " (math round [det] 5) " from the U diagonal, matrix det gives " (math round (matrix det [M]) 5) "; output x = " [x] "; multipliers " [mult]
Result
mentdb
det(I-A) = 0.36287 from the U diagonal, matrix det gives 0.36287; output x = [280.4684808818464,763.4860489149157,840.9920771615571]; multipliers [2.1357216672407855,2.4733034791594903,1.9703754736479504]
Go further
  • Check the Hawkins-Simon condition: all leading minors of I - A positive
  • Simulate a 10 % rise of the services demand
executed on MentDB

Discrete newsvendor from an empirical demand

A bakery knows the daily sales of a fresh cake over 20 days. The script builds the frequency table, the critical ratio Cu/(Cu + Co) and picks the smallest quantity whose cumulative frequency reaches it.

Pipeline
01 · tableFrequency tablestat frequency→02 · ratioCritical ratiomath calc→03 · pickFirst cumulative frequency >= ratiojson select→04 · modeMost frequent demandstat mode
Script
admin
#1. Daily sales of a fresh cake over 20 days: cost 4 EUR, price 10 EUR, unsold sold off at 1 EUR;
-> "[sales]" "[18, 22, 20, 25, 20, 19, 23, 22, 20, 26, 21, 18, 22, 24, 20, 22, 19, 23, 21, 22]";
#2. Empirical distribution of the demand;
json load "rx_p4b_freq" (stat frequency [sales]);
#3. Critical ratio Cu / Cu + Co: under-stock costs 10-4, over-stock costs 4-1;
-> "[cr]" (math calc "(10-4)/((10-4)+(4-1))");
#4. Order the smallest quantity whose cumulative frequency reaches the critical ratio;
-> "[q]" ""; -> "[i]" 0;
while (is empty [q]) {
	if (>= (json select "rx_p4b_freq" (concat "/[" [i] "]/cumulative_frequency")) [cr]) { -> "[q]" (json select "rx_p4b_freq" (concat "/[" [i] "]/value")); };
	++ "[i]";
};
-> "[nv]" (json count "rx_p4b_freq" /);
json unload "rx_p4b_freq";
concat [nv] " distinct demand values, most frequent " (stat mode [sales]) ", mean " (stat mean [sales]) "; critical ratio " (math round [cr] 4) ": bake " [q] " cakes a day"
Result
mentdb
9 distinct demand values, most frequent [22], mean 21.35; critical ratio 0.6667: bake 22 cakes a day
Go further
  • Compute the expected profit of each quantity from the table
  • Compare with the normal approximation stat dist_inv

Also in Sciences

MA
Mathematician
Calculus, linear algebra, number theory and complex numbers: exact and numerical in the same script.
10 recipes
MT
Math teacher & student
Equations, tables of values, fractions and probabilities: exact answers and their decimal check.
10 recipes
PH
Physicist
Mechanics with units, constants and differential equations: trajectories, pendulums, oscillators.
10 recipes
OP
Optics & photonics
Refraction, lenses, photons and diffraction with SI units and CODATA constants.
10 recipes
TH
Thermodynamics & heat transfer
Conduction, heat pumps, ideal gases and cooling laws with units, integrals and ODEs.
10 recipes
NR
Nuclear & radiation physicist
Decay, half-lives, generators and shielding: dates, fits and solvers for radiation work.
10 recipes
AQ
Atomic & quantum physicist
Spectral lines, matter waves, blackbody radiation and CODATA constants with their uncertainties.
10 recipes
AE
Acoustics engineer
Decibels, Doppler, spectra and reverberation: the daily formulas of noise and room acoustics.
10 recipes
CN
Cryptographer & number theorist
RSA and Diffie-Hellman by hand, primes, modular arithmetic and exact big integers.
10 recipes
ME
Metrologist
Uncertainty budgets following the GUM, exact decimals and calibration curves.
10 recipes
MP
Mathematical physicist
Special functions at work: Bessel modes, Hermite and Laguerre states, gamma and zeta, with the CAS.
10 recipes
CM
Crystallographer & materials physicist
Miller indices, metric tensors, elastic tensors and structure factors.
10 recipes