Skip to content

Instantly share code, notes, and snippets.

@bqpd

bqpd/Optimizing Maintenance Schedules.ipynb

Last active December 22, 2015 16:58
Show Gist options
  • Save bqpd/6502461 to your computer and use it in GitHub Desktop.
Save bqpd/6502461 to your computer and use it in GitHub Desktop.
Download ZIP
see it on the iPython Notebook Viewer at http://nbviewer.ipython.org/6502461
Display the source blob
Display the rendered blob
Raw
Optimizing Maintenance Schedules.ipynb
{
"metadata": {
"name": ""
},
"nbformat": 3,
"nbformat_minor": 0,
"worksheets": [
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Optimizing Maintenance Schedules\n",
"--------------------------------\n",
"\n",
"_Updated Sept 9th, 2013 by Ned Burnell. Find an up to date version at bit.ly/maintop_\n",
"\n",
"This script minimizes the maintenance costs of a system which needs to ensure a certain system reliability in the face of multiple independently-failing parts. To do so it phrases the problem as one of convex optimization, and solves it with the [cvxopt](cvxopt.org) package.\n",
"\n",
"\n",
"**CVXopt Installation Instructions**\n",
"\n",
" 1. Download the [gzipped tar file](https://github.com/cvxopt/cvxopt/archive/1.1.6.tar.gz)\n",
" 2. Decompress it\n",
" 3. Install it\n",
" \n",
"In linux or Mac OS X, you can do all this by running the following in a terminal:\n",
"\n",
" wget https://github.com/cvxopt/cvxopt/archive/1.1.6.tar.gz\n",
" tar zxvf 1.1.6.tar.gz\n",
" cd cvxopt-1.1.6\n",
" sudo python setup.py install\n",
"\n",
"\n",
"**Preamble**\n",
"\n",
"The first part of this script just imports the required functions and modules:"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"from cvxopt import matrix\n",
"from cvxopt.solvers import cpl, options\n",
"from math import log, log10"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 1
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"System Properties\n",
"-----------------\n",
"\n",
"$L$ is the desired system lifetime.\n",
"\n",
"$R_{min}$ is the minimum allowable system reliability.\n",
"\n",
"$N$ is the number of parts to be modeled in the system.\n",
"\n",
"**Part properties**\n",
"\n",
"$c =\\left(\\begin{array}{ccc}\n",
"c_1 \\\\\n",
"\\vdots \\\\\n",
"c_N \\end{array}\\right)$ is the cost to replace each of those parts, in arbitrary units.\n",
"\n",
"Each part is currently assumed to have a polynomial failure rate, corresponding to a Weibull distribution, but any distribution that meets certain criteria (see the Feasibility Note in the Log-Reliability section) could be optimized in exactly the same way.\n",
"\n",
"Constant failure rates (exponential distributions) are obviously not worth replacing before they fail, and any such parts will effectively just change $R_{min}$.\n",
"\n",
"*Weibull distribution properties*\n",
"\n",
"$L_{10} = \\left(\\begin{array}{ccc}\\vdots\\end{array}\\right)$ is the cost to replace each of those parts, in arbitrary units.\n",
"\n",
"$\\Lambda = 4.828 L_{10}$ is the scale factor for the ISO-281 Weibull distribution.\n",
"\n",
"$k = \\left(\\begin{array}{ccc}\\vdots\\end{array}\\right)$ is the shape factor for each Weibull distribution."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"## System Properties\n",
"L = 0.1 # desired system lifetime\n",
"Rmin = 0.997 # minimum system reliability\n",
"N = 3 # number of parts in the system\n",
"\n",
"## Part properties\n",
"c = matrix([1., 2., 2.]) # cost to replace each part\n",
"L10 = matrix([1, 3, 1]) # L10 lifetimes, in same units as L\n",
"lam = 4.4828*L10 # scale factors of their Weibull distributions\n",
"k = matrix([1.5, 1.5, 2]) # shape factors (1.5 is what's used by ISO-281)\n",
"\n",
"## Check that the dimensions are correct\n",
"for vector in (c, L10, lam, k): assert(vector.size == (N, 1))\n",
" # IF THE ASSERT FAILS, CHECK WHETHER THE VECTORS ARE ALL THE RIGHT SIZE"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 2
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Log-Reliability\n",
"---------------\n",
"\n",
"Before we can give the problem to the solver, we need a function which can return the log-reliability of the system given a particular maintenance schedule, as well as its first and second derivatives.\n",
"\n",
"$x_i$ is the frequency at which part $i$ is replaced, in the inverse units of $L$. (e.g., if $L = 1$, part $i$ will be replaced $x_i$ times)\n",
"\n",
"$p_i$ is the probability distribution of part $i$ failing.\n",
"\n",
"$r_i = 1 - \\int_0^\\frac{1}{x_i}{p_i(t)dt}$ is the chance that part $i$ will not fail before being replaced.\n",
"\n",
"$R = \\prod_{i=1}^N r_i^{L x_i}$ is the system's reliability (the probability that it will not have failed by time $L$)\n",
"\n",
"$Q = \\frac{log(R)}{L} = \\sum_{i=1}^N x_i\\log{r_i}$ is the lifetime-adjusted log-reliability.\n",
"\n",
"$q_i = x_i\\log{r_i}$ is the part-specific lifetime-adjusted log-reliability.\n",
"\n",
"$Q_{min} = \\frac{\\log(R_{min})}{L}$ is the lifetime-adjusted log-reliability constraint.\n",
"\n",
"_Weibull distribution specific_\n",
"\n",
"$q_i^{W} = -\\lambda_i^{-k_i} x_i^{-k_i+1}$ is $q_i$ for a Weibull distribution of shape $k_i$ and scale $\\lambda_i$\n",
"\n",
"*Feasibility Note*\n",
"\n",
"For the problem to be feasible, the Hessian needs to be negative semi-definite for all possible values of $x$: for this problem that's equivalent to saying that for all $i$, $\\frac{d^2 q_i}{dx_i^2} \\leq 0$.\n",
"\n",
"\n",
"**Log-reliability function**\n",
"\n",
"The code below generates the function $F(x) = Q_{min} - Q(x)$ from the parameters of $\\Lambda, k, R_{min}$ and $L$.\n",
"\n",
"Note that the function will also return derivatives and the Hessian of $F$ when the solver asks for them, and multiply the Hessian by a solver-provided scale factor.\n",
"\n",
"When the function is first generated it prints out its initial guess for the solution $x^*$: a useful hack if the optimizer isn't working is to tweak constants in the functions that generate that guess until they're close to the optimizer's unfinished guess."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"def F_gen(N, lam, k, Rmin, L):\n",
" \"\"\"Generates function F(x) = Q_min - Q(x)\"\"\"\n",
" Qmin = log(Rmin)/L\n",
" \n",
" # Scaling parameters\n",
" magQ = log10(-Qmin)+0.6\n",
" s = 10**(-4*magQ)\n",
" x0 = matrix(10**(-2*magQ),(N,1))\n",
" print \"x0 is\", x0.T\n",
" \n",
" # Log-reliability function to return\n",
" def log_reliability_function(x=None, z=None):\n",
" if x is None: \n",
" return 1, x0\n",
" else:\n",
" # all elements of x must be nonnegative\n",
" for e in x:\n",
" if e <= 0: return None\n",
" \n",
" f = Qmin\n",
" Df = matrix(0.,(1,N))\n",
" H = matrix(0.,(N, N))\n",
" \n",
" for i in range(N):\n",
" f = f + lam[i]**-k[i] * x[i]**(-k[i]+1)\n",
" Df[i] = lam[i]**-k[i] * x[i]**(-k[i]) * (-k[i]+1)\n",
" H[i,i] = lam[i]**-k[i] * x[i]**(-k[i]-1) * (-k[i]+1) * -k[i]\n",
" \n",
" if z is not None: return s*f, s*Df, s*z*H\n",
" else: return s*f, s*Df\n",
" \n",
" return log_reliability_function\n",
"\n",
"F = F_gen(N, lam, k, Rmin, L)"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"x0 is [ 6.99e+01 6.99e+01 6.99e+01]\n",
"\n"
]
}
],
"prompt_number": 3
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Using a Convex-Optimization Solver\n",
"----------------------------------\n",
"\n",
"Since we know the cost ($c$) to replace each part and require a minimum system reliability ($R_{min}$), minimize the expected maintenance cost of the system ($c^Tx$) by changing how frequently each part is replaced ($x$), making sure that the resulting maintenance schedule makes the system at least as reliable as it needs to be ($R(x) \\geq R_{min}$). In the language of optimization this looks like:\n",
"\n",
"$\\begin{aligned}\n",
"& \\underset{x}{\\text{minimize}}\n",
"& & c^Tx \\\\\n",
"& \\text{subject to}\n",
"& & R(x) \\geq R_{min} \\\\\n",
"&&& x \\gt 0.\n",
"\\end{aligned}$\n",
"\n",
"The constraint $ R(x) \\geq R_{min} $ can, by taking the log of each side and dividing by $L$, be rewritten as $ Q(x) \\geq Q_{min} $, which can be rearranged as $ Q_{min} - Q(x) = F(x) \\leq 0 $ to fit the form we need for the <tt>cpl</tt> solver of <tt>cvxopt</tt>.\n",
"\n",
"Run <tt>help(cpl)</tt> at a python prompt to learn more about the input requirements of <tt>cpl</tt>."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"options['maxiters'] = 20 # maximum iterations to allow\n",
"options['feastol'] = 1e-7 # minimum accuracy of the final solution\n",
"solution = cpl(c,F)"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
" pcost dcost gap pres dres\n",
" 0: 3.4948e+02 2.7959e+02 1e+00 1e+00 1e+00\n",
" 1: 1.2649e+02 1.8462e+02 9e-01 7e-01 2e+01\n",
" 2: 8.4704e+01 1.4342e+02 5e-01 6e-01 2e+01\n",
" 3: 8.0134e+01 1.0430e+02 6e-03 2e-01 7e+00\n",
" 4: 7.5661e+01 8.2902e+01 9e-04 9e-02 3e+00\n",
" 5: 7.0135e+01 7.3932e+01 1e-04 6e-02 9e-01\n",
" 6: 6.9527e+01 7.1041e+01 1e-05 3e-02 3e-01\n",
" 7: 7.0282e+01 7.0589e+01 4e-07 5e-03 4e-02\n",
" 8: 7.0560e+01 7.0574e+01 4e-09 2e-04 2e-03\n",
" 9: 7.0573e+01 7.0574e+01 4e-11 3e-06 2e-05\n",
"10: 7.0574e+01 7.0574e+01 4e-13 3e-08 2e-07"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n",
"11: 7.0574e+01 7.0574e+01 4e-15 3e-10 2e-09\n",
"Optimal solution found.\n"
]
}
],
"prompt_number": 4
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"solution"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 5,
"text": [
"{'dual infeasibility': 1.970110761655239e-09,\n",
" 'dual objective': 70.57355075192288,\n",
" 'dual slack': 0.825136701012636,\n",
" 'gap': 4.402332001590911e-15,\n",
" 'primal infeasibility': 2.8395488425850437e-10,\n",
" 'primal objective': 70.57355073578202,\n",
" 'primal slack': 5.335275956321199e-15,\n",
" 'relative gap': 6.237934686134469e-17,\n",
" 'sl': <0x1 matrix, tc='d'>,\n",
" 'snl': <1x1 matrix, tc='d'>,\n",
" 'status': 'optimal',\n",
" 'x': <3x1 matrix, tc='d'>,\n",
" 'y': <0x1 matrix, tc='d'>,\n",
" 'zl': <0x1 matrix, tc='d'>,\n",
" 'znl': <1x1 matrix, tc='d'>}"
]
}
],
"prompt_number": 5
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"# Print out the final result\n",
"print \"The minimal cost of\", solution['dual objective']\n",
"print \" is achieved at a maintenance frequency of \\n \", solution['x'].T"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"The minimal cost of 70.5735507519\n",
" is achieved at a maintenance frequency of \n",
" [ 3.56e+01 7.47e+00 1.00e+01]\n",
"\n"
]
}
],
"prompt_number": 6
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"TODO\n",
"----\n",
"\n",
" 1. by running many trials, determine where the solver is stable, and try to make it consistently stable for different shape parameters / number of parts / desired reliabilities\n",
" 2. make a script to automatically determine the cost of replacing one part with another, and use it to analyze various servo options\n",
" 3. make sensitivity plots and try alternate solvers, like Newton's method"
]
}
],
"metadata": {}
}
]
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment