In [1]:

    

%pylab inline
from gpkit import GP, VectorVariable, Variable
import gpkit.interactive
gpkit.interactive.init_printing()


    

    




Populating the interactive namespace from numpy and matplotlib

In [2]:

    

L = Variable("L", 0.1, "years", "system lifetime")

Rmin = 0.997  # minimum system reliability
mlogRmin = Variable("\\mathrm{log}(\\tfrac{1}{R_{min}})",
                    log(1/Rmin), "-", "log of inverse of system reliability")

n = 3  # number of parts

c = VectorVariable(n, "c", [1, 2, 2], "-", "cost to replace each part, in USD")

q = gpkit.VectorVariable(n, "q", "1/years", "log of inverse of part reliabilities")
L10 = VectorVariable(n, "L_{10}", [1, 3, 1], "years", "L10 lifetimes")
k = array([1.5, 1.5, 2])  # part shape factors; 1.5 is what's used by ISO-281
lam = 4.4828*L10  # Weibull scale factors

x = gpkit.VectorVariable(n, "x", "1/years", "maintenance rate")


    

In [3]:

    

gp = gpkit.GP(L*(c*x).sum(),
              [q == x*(lam*x)**-k,
               q.sum() <= mlogRmin/L,
               L*x >= 1])
gp


    

    
Out[3]:





$$\begin{array}[ll]
\text{}
\text{minimize}
    & {c}_{0} {x}_{0} L + {x}_{1} {c}_{1} L + {x}_{2} {c}_{2} L\mathrm{\left[ - \right]} \\
\text{subject to}
    & {q}_{0} = \frac{0.1054}{{L_{10}}_{0}^{1.5} {x}_{0}^{0.5}} \\
    & {q}_{1} = \frac{0.1054}{{x}_{1}^{0.5} {L_{10}}_{1}^{1.5}} \\
    & {q}_{2} = \frac{0.04976}{{x}_{2} {L_{10}}_{2}^{2}} \\
    & \frac{\mathrm{log}(\tfrac{1}{R_{min}})}{L} \geq {q}_{0} + {q}_{1} + {q}_{2} \\
    & 1 \leq {x}_{0} L \\
    & 1 \leq {x}_{1} L \\
    & 1 \leq {x}_{2} L \\
\text{substituting}
    & L = 0.1 \\
    & \mathrm{log}(\tfrac{1}{R_{min}}) = 0.003005 \\
    & {L_{10}}_{0} = 1 \\
    & {L_{10}}_{1} = 3 \\
    & {L_{10}}_{2} = 1 \\
    & {c}_{0} = 1 \\
    & {c}_{1} = 2 \\
    & {c}_{2} = 2 \\
\end{array}$$

In [4]:

    

sol = gp.solve()
print 
print "The optimal schedule, over the lifetime of %.2g years" % sol(L)
print "    is to buy each part [ %s ] times" % "  ".join(["%.2g" % xi for xi in sol(L*x)])
print "    at a total cost of $%.3g" % sol["cost"]


    

    




Using solver 'cvxopt'
Solving for 6 variables.
Solving took 0.0813 seconds.

The optimal schedule, over the lifetime of 0.1 years
    is to buy each part [ 3.2  1  1 ] times
    at a total cost of $7.19

In [5]:

    

sol["local_model"].flatten()[0]


    

    
Out[5]:





$$0.3675\frac{{c}_{0}^{0.44} {c}_{1}^{0.28} {c}_{2}^{0.28} L^{2.3}}{\mathrm{log}(\tfrac{1}{R_{min}})^{1.4} {L_{10}}_{0}^{1.3} {L_{10}}_{1}^{0.46} {L_{10}}_{2}^{0.47}}$$

In [6]:

    

from IPython import utils
from IPython.core.display import HTML
import os
def css_styling():
    """Load default custom.css file from ipython profile"""
    base = utils.path.get_ipython_dir()
    styles = ("<style>\n%s\n</style>" % 
              open(os.path.join(base,'profile_default/static/custom/custom.css'),'r').read())
    return HTML(styles)
css_styling()


    

    
Out[6]: