denis-bz/0-Lower-risk-stochastic-programming-Dakota-problem.md

## 0-Lower-risk-stochastic-programming-Dakota-problem.md

      
    Raw
  

              0-Lower-risk-stochastic-programming-Dakota-problem.md
            
          
    Easy lower-risk stochastic programming: the Dakota problem

Keywords: stochastic programming, risk, linear programming, python
Users of stochastic programming sometimes look only at "expected" payoff and ignore risk.
For example, a well-known tutorial problem, Dakota furniture from Higle 2005,
gives a max-expected profit $ 1730,
but with 30 % chance of $ 650 loss, 70 % $2750 profit.
That looks risky to me (an engineer, not a businessman).
The profit - risk tradeoff


It's easy to run SP for less risk:
constrain profit to be >= for example $ 0, no loss.
4 runs of Dakota LP with constraints None, $ 0, $ 500, and $ 1000
give this table of profit when demand is low, medium or high:
          low  medium  high demand
          30   40   30 % probability
------------------------------------------------------
profit $ [-650 2750 2750]  average $ 1730  cost $ 6450  -- Higle
profit $ [   0 2415 2415]  average $ 1691  cost $ 5800  -- constrain profit >= $ 0
profit $ [ 500 1902 1902]  average $ 1481  cost $ 4581  -- constrain profit >= $ 500
profit $ [1000 1342 1342]  average $ 1240  cost $ 3234  -- constrain profit >= $ 1000
------------------------------------------------------

Summary


don't optimize expected payoff alone, look at risk too
Dakota is only one (1) textbook problem; I'd appreciate links to more.

The Dakota furniture problem

Dakota furniture makes desks, tables and chairs on a monthly schedule, as follows:

on the 1st of each month, buy lumber etc. (3 resources) for the whole month
then get a forecast whether demand for the month will be low, medium or high.
Based on that, decide how many desks, how many tables, and how many chairs
to make and sell over the month.
Desks sell for $60, tables $40, chairs $10.
There are resource limits: Dakota bought only so much lumber etc. on the 1st.

(Of course, this model is unrealistically simple, but.)
The "recourse" stochastic programming solution on Higle p. 36 is:

buy lumber etc. on each 1st for $1300 + $540 + $325, then make
     Desks Tables Chairs  Profit
------------------------------
low:    50   20  200     $ -650 (loss)
medium: 80  110    0     $ 2750
high:   80  110    0     $ 2750

(The general problem of stochastic programming is to decide what and how much to buy now,
to meet uncertain demand over the next say 7 days or 30 days or 12 months.
Examples: chain store inventory, cars, buildings.
There are many approaches to such a range of problems,
with subcultures who think differently: business, finance, optimization ...)
Program outline

dakota_abc.py: build LP matrices A b c in numpy or scipy.sparse format

dakota.py: solve with
scipy linprog interior-point

(Run python dakota.py T=2 to see the structure of A.)
The parameter T= builds larger matrices, 3T+1 x 3T+3, for T time periods.
But there's no coupling between time periods, so solutions are independent of T.
Notes

There are many papers on Benders etc. decomposition schemes
that claim to speed up big stochastic-programming LPs.
How do these compare with sparse interior-point linprog,
on other opensource problems on the web ?
See also

Higle 2005, "Stochastic Programming: Optimization When Uncertainty Matters", 24 pages,
pdf

Wikipedia Stochastic programming

github.com/topics/stochastic-optimization

splib test cases -- rather a tangle
Comments welcome, test cases most welcome.

cheers

-- denis 29 May 2019
_{Theory and practice are closer in theory than they are in practice.}

  
## dakota.py
#!/usr/bin/env python
""" expand the Dakota stochastic programming problem from Higle 2005 p. 35 """
# see https://en.wikipedia.org/wiki/Stochastic_programming#Deterministic_equivalent_of_a_stochastic_problem
# linprog == Higle low medium high  for T > 1 too

from __future__ import division, print_function
import resource
import sys
import numpy as np
from numpy import r_, c_, inf
from scipy.optimize import linprog  # $scopt/_linprog.py linprog.doc
    # https://docs.scipy.org/doc/scipy/reference/optimize.linprog-interior-point.html

from dakota_abc import dakota_abc
import znumpyutil as nu

np.set_printoptions( threshold=20, edgeitems=11, linewidth=140,
        formatter = dict( float=nu.numstr ))
        # formatter = dict( float = lambda x: "%.2g" % x ))  # float arrays %.2g
print("\n" + 80 * "=")

#...........................................................................
    # Dakota problem parameters --
T = 1  # > 1: stack A (3T, 3 + 3T)  T=100: A (900, 903)
probs = r_[ 3., 4, 3 ]
ccost = r_[ 2., 4, 5.2 ]
csell = r_[ 60., 40, 10 ]

a = np.array([  # a . [sell D T C] <= resources x
    # D T C
    [8., 6, 1],     # <= Lumber aka x_b
    [4, 2, 1.5],    # <= Finishing x_f
    [2, 1.5, 0.5] ]) # <= Carpentry x_c

demand = r_[
    # D T C
     50., 20, 200,  # low
    150, 110, 225,  # medium
    250, 250, 500 ] # high

save = 0

    # linprog --
sparse = False  # False: print A
disp = True  # patch $scipy/optimize/_linprog_ip.py clock, alpha NaN not "-"
maxiter = 1e3
tol = 1e-6  # default 1e-8

minprofit = - inf  # T=1 dense

    # to change these params, run this.py  a=1  b=None  'c = expr' ...
    # in sh or ipython --
for arg in sys.argv[1:]:
    exec( arg )

probs = np.asarray( probs, dtype=float )
probs /= probs.sum()

params = "%s  T %d  probs %s  ccost %s  csell %s  tol %g  sparse %s  minprofit %g \n" % (
        sys.argv[0], T, probs, ccost, csell, tol, sparse, minprofit )
print("#", params)


def dakota_print( x, c, ccost=ccost, csell=csell, probs=probs ):
    """ unpack, print cost profit x c*x """
    ncost = len(ccost)
    nsell = len(csell)
    cost = x[:ncost] .dot( ccost )
    nsell = x[ncost:].reshape( -1, nsell )
    profit = (nsell.dot( csell ) - cost) .reshape( -1, len(csell) )
    av = profit[0].dot( probs )

    print("\n" + 80 * "-")
    print("cost: $ %.0f = %s * $ %s  lumber, finishing and carpentry " % (
            cost, x[:ncost], ccost ))
    print("         low medium high %s %%" % (probs * 100))
    print("profit $ %s  average $ %.0f  cost $ %.0f " % (
            nu.ints( profit[0] ), av, cost ))
    print("number of desks tables chairs sold: \n", nu.ints( nsell.T ))
    print(80 * "-")

    cmax = - c
    print("cmax:", cmax)
    print("x:   ", x)
    print("c*x: ", cmax * x)


#...............................................................................
    # expand the Dakota stochastic programming problem from Higle 2005 p. 35
    # -> A 3T, 3T+3 for linprog
A, b, c = dakota_abc( a=a, ccost=ccost, csell=csell, probs=probs, T=T,
        sparse=sparse, minprofit=minprofit )
hibounds = r_[ len(ccost) * [np.inf], T * tuple( demand )]
bounds = c_[ np.zeros_like( hibounds ), hibounds ]
print("bounds: lumber finishing carpentry   D T C low, medium, high demand")
print(hibounds)
options = dict( disp=disp, maxiter=int(maxiter), sparse=sparse, tol=tol )
print("\nscipy linprog interior-point --")

#...............................................................................
res = linprog( c, A_ub=A, b_ub=b, method="interior-point", bounds=bounds,
        options=options )
x = res.x
if disp is False:
    if res.status != 0:
        print(res.message)
    print("linprog: -f %g " % (- res.fun))

u = resource.getrusage( resource.RUSAGE_SELF )  # man 2 getrusage
print("max memory: %.0f M" % (u.ru_maxrss / 1e6))

dakota_print( x, c )

# b_Ax = b - A.dot( x )  # >= 0 ?
# print "b - Ax:", b_Ax

if save:
    out = "%d-dakota.npz" % T
    print("\nsaving to", out)
    np.savez( out, A=A, b=b, c=c, x=x, T=T, params=params )

## dakota_abc.py
#!/usr/bin/env python
""" expand the Dakota stochastic programming problem from Higle 2005 p. 35 """

from __future__ import division, print_function
import numpy as np
from numpy import r_, c_, inf

from znumpyutil import blank0

#...............................................................................
def dakota_abc( a, ccost, csell, probs, T=1, sparse=True, minprofit=-inf ):
    """ expand the Dakota stochastic programming problem from Higle 2005 p. 35
        -> A 3T, 3T+3  for linprog( c, A_ub=A, b_ub=b )
    """
    from scipy.linalg import block_diag  # (*args)
    if sparse:
        from scipy import sparse as sp  # eye hstack vstack
        # from scipy.sparse import block_diag  # (tup)
    else:
        sp = np

    cprobsell = np.outer( probs, csell ).reshape(-1)  # [p0 csell, p1 csell ...]
    c = r_[ - T * ccost, T * tuple( cprobsell )]  # maximize
    Tp = T * len(probs)
    I = - sp.eye( a.shape[0] )
    minusIs = sp.vstack(( Tp * [I] ))
    blocks = block_diag( *( Tp * [a] ))  # a, a, a  sparse ?
    A = sp.hstack(( minusIs, blocks ))

    if minprofit > -inf:  # add constraint - ccost . x + csell . y0 >= minprofit
        assert T == 1 and not sparse
        row = np.zeros_like( A[0] )
        ncost = len(ccost)
        nsell = len(csell)
        row[:ncost] = ccost
        row[ncost : ncost+nsell] = - csell
        A = np.vstack(( A, row ))

    if sparse:
        A = A.tocsc()  # ?
        nbytes = A.data.nbytes + A.indices.nbytes + A.indptr.nbytes
        nnz = A.nnz
    else:
        nbytes = A.nbytes
        nnz = np.count_nonzero(A)

    print("dakota_Abc T=%d: A %s  %.3g bytes  %d nnz" % (
            T, A.shape, nbytes, nnz ))
    print("""
buy lumber finishing carpentry,  make and sell desks tables chairs
  L  F  C   D  T  C low demand
                     D  T  C medium demand
                              D  T  C high demand
""")
    if not sparse:
        print(blank0(A))
    else:
        print(blank0( A[:4].A ), "... sparse")
    print("c:", c)
    b = np.zeros( A.shape[0] )
    if minprofit > -inf:
        b[-1] = - minprofit
    return A, b, - c


## minprofit1000-dakota.log
# from: dakota.py T=1 minprofit=1000 sparse=0
# 24 May 2019 10:02 gmt  ~bz/py/opt/lp/linprog  Denis-iMac 10.10.5
versions: numpy 1.16.3  scipy 1.2.1  python 2.7.15  mac 10.10.5


================================================================================
# dakota.py  T 1  probs [0.3 0.4 0.3]  ccost [2 4 5.2]  csell [60 40 10]  tol 1e-06  sparse 0  minprofit 1000

dakota_Abc T=1: A (10, 12)  960 bytes  42 nnz

buy lumber finishing carpentry,  make and sell desks tables chairs
  L  F  C   D  T  C low demand
                     D  T  C medium demand
                              D  T  C high demand

[[ -1           8   6   1                        ]
 [     -1       4   2   1                        ]
 [         -1   2   1                            ]
 [ -1                       8   6   1            ]
 [     -1                   4   2   1            ]
 [         -1               2   1                ]
 [ -1                                   8   6   1]
 [     -1                               4   2   1]
 [         -1                           2   1    ]

 [  2   4   5  -60 -40 -10                       ]]  -- cost - profit >= 1000

c: [-2 -4 -5.2  18 12 3    24 16 4     18 12 3]
                30         40          30 % * $ 60 40 10

bounds: lumber finishing carpentry   D T C low, medium, high demand
[inf inf inf  50 20 200  150 110 225  250 250 500]

scipy linprog interior-point --
Sec  Objective    Primal, Dual Feasibility   Gap   Step   Path Parameter
  0: -98.8        1          1          1          nan        1
  0: -96.6318     0.432      0.432      0.432      0.578      0.432
  0: -108.421     0.113      0.113      0.113      0.761      0.113
  0: -154.069     0.0292     0.0292     0.0292     0.75       0.0292
  0: -486.968     0.00597    0.00597    0.00597    0.813      0.00597
  0: -991.449     0.0013     0.0013     0.0013     0.803      0.0013
  0: -1223.89     8.92e-05   8.92e-05   8.92e-05   0.963      8.92e-05
  0: -1236.61     1.49e-05   1.49e-05   1.49e-05   0.846      1.49e-05
  0: -1239.65     3.33e-08   3.33e-08   3.33e-08   1          3.33e-08
Optimization terminated successfully.
         Current function value: -1239.653881
         Iterations: 8
max memory: 39 M

--------------------------------------------------------------------------------
cost: $ 3234 = [610 305 153] * $ [2 4 5.2]  lumber, finishing and carpentry
         low medium high [30 40 30] %
profit $ [1000 1342 1342]  average $ 1240  cost $ 3234
number of desks tables chairs sold:
[[50 76 76]
 [20  0  0]
 [43  0  0]]
--------------------------------------------------------------------------------
cmax: [-2 -4 -5.2 18 12 3 24 16 4 18 12 3]
x:    [610 305 153  50 20 43  76 0.00018 0.00061  76 0.00024 0.00062]
c*x:  [-1220 -1220 -793 900 240 130 1830 0.0028 0.0025 1373 0.0029 0.0019]

        0.48 real         0.86 user         0.28 sys
	#!/usr/bin/env python
	""" expand the Dakota stochastic programming problem from Higle 2005 p. 35 """
	# see https://en.wikipedia.org/wiki/Stochastic_programming#Deterministic_equivalent_of_a_stochastic_problem
	# linprog == Higle low medium high for T > 1 too

	from __future__ import division, print_function
	import resource
	import sys
	import numpy as np
	from numpy import r_, c_, inf
	from scipy.optimize import linprog # $scopt/_linprog.py linprog.doc
	# https://docs.scipy.org/doc/scipy/reference/optimize.linprog-interior-point.html

	from dakota_abc import dakota_abc
	import znumpyutil as nu

	np.set_printoptions( threshold=20, edgeitems=11, linewidth=140,
	formatter = dict( float=nu.numstr ))
	# formatter = dict( float = lambda x: "%.2g" % x )) # float arrays %.2g
	print("\n" + 80 * "=")

	#...........................................................................
	# Dakota problem parameters --
	T = 1 # > 1: stack A (3T, 3 + 3T) T=100: A (900, 903)
	probs = r_[ 3., 4, 3 ]
	ccost = r_[ 2., 4, 5.2 ]
	csell = r_[ 60., 40, 10 ]

	a = np.array([ # a . [sell D T C] <= resources x
	# D T C
	[8., 6, 1], # <= Lumber aka x_b
	[4, 2, 1.5], # <= Finishing x_f
	[2, 1.5, 0.5] ]) # <= Carpentry x_c

	demand = r_[
	# D T C
	50., 20, 200, # low
	150, 110, 225, # medium
	250, 250, 500 ] # high

	save = 0

	# linprog --
	sparse = False # False: print A
	disp = True # patch $scipy/optimize/_linprog_ip.py clock, alpha NaN not "-"
	maxiter = 1e3
	tol = 1e-6 # default 1e-8

	minprofit = - inf # T=1 dense

	# to change these params, run this.py a=1 b=None 'c = expr' ...
	# in sh or ipython --
	for arg in sys.argv[1:]:
	exec( arg )

	probs = np.asarray( probs, dtype=float )
	probs /= probs.sum()

	params = "%s T %d probs %s ccost %s csell %s tol %g sparse %s minprofit %g \n" % (
	sys.argv[0], T, probs, ccost, csell, tol, sparse, minprofit )
	print("#", params)


	def dakota_print( x, c, ccost=ccost, csell=csell, probs=probs ):
	""" unpack, print cost profit x c*x """
	ncost = len(ccost)
	nsell = len(csell)
	cost = x[:ncost] .dot( ccost )
	nsell = x[ncost:].reshape( -1, nsell )
	profit = (nsell.dot( csell ) - cost) .reshape( -1, len(csell) )
	av = profit[0].dot( probs )

	print("\n" + 80 * "-")
	print("cost: $ %.0f = %s * $ %s lumber, finishing and carpentry " % (
	cost, x[:ncost], ccost ))
	print(" low medium high %s %%" % (probs * 100))
	print("profit $ %s average $ %.0f cost $ %.0f " % (
	nu.ints( profit[0] ), av, cost ))
	print("number of desks tables chairs sold: \n", nu.ints( nsell.T ))
	print(80 * "-")

	cmax = - c
	print("cmax:", cmax)
	print("x: ", x)
	print("cx: ", cmax x)


	#...............................................................................
	# expand the Dakota stochastic programming problem from Higle 2005 p. 35
	# -> A 3T, 3T+3 for linprog
	A, b, c = dakota_abc( a=a, ccost=ccost, csell=csell, probs=probs, T=T,
	sparse=sparse, minprofit=minprofit )
	hibounds = r_[ len(ccost) * [np.inf], T * tuple( demand )]
	bounds = c_[ np.zeros_like( hibounds ), hibounds ]
	print("bounds: lumber finishing carpentry D T C low, medium, high demand")
	print(hibounds)
	options = dict( disp=disp, maxiter=int(maxiter), sparse=sparse, tol=tol )
	print("\nscipy linprog interior-point --")

	#...............................................................................
	res = linprog( c, A_ub=A, b_ub=b, method="interior-point", bounds=bounds,
	options=options )
	x = res.x
	if disp is False:
	if res.status != 0:
	print(res.message)
	print("linprog: -f %g " % (- res.fun))

	u = resource.getrusage( resource.RUSAGE_SELF ) # man 2 getrusage
	print("max memory: %.0f M" % (u.ru_maxrss / 1e6))

	dakota_print( x, c )

	# b_Ax = b - A.dot( x ) # >= 0 ?
	# print "b - Ax:", b_Ax

	if save:
	out = "%d-dakota.npz" % T
	print("\nsaving to", out)
	np.savez( out, A=A, b=b, c=c, x=x, T=T, params=params )
	# from: dakota.py T=1 minprofit=1000 sparse=0
	# 24 May 2019 10:02 gmt ~bz/py/opt/lp/linprog Denis-iMac 10.10.5
	versions: numpy 1.16.3 scipy 1.2.1 python 2.7.15 mac 10.10.5


	================================================================================
	# dakota.py T 1 probs [0.3 0.4 0.3] ccost [2 4 5.2] csell [60 40 10] tol 1e-06 sparse 0 minprofit 1000

	dakota_Abc T=1: A (10, 12) 960 bytes 42 nnz

	buy lumber finishing carpentry, make and sell desks tables chairs
	L F C D T C low demand
	D T C medium demand
	D T C high demand

	[[ -1 8 6 1 ]
	[ -1 4 2 1 ]
	[ -1 2 1 ]
	[ -1 8 6 1 ]
	[ -1 4 2 1 ]
	[ -1 2 1 ]
	[ -1 8 6 1]
	[ -1 4 2 1]
	[ -1 2 1 ]

	[ 2 4 5 -60 -40 -10 ]] -- cost - profit >= 1000

	c: [-2 -4 -5.2 18 12 3 24 16 4 18 12 3]
	30 40 30 % * $ 60 40 10

	bounds: lumber finishing carpentry D T C low, medium, high demand
	[inf inf inf 50 20 200 150 110 225 250 250 500]

	scipy linprog interior-point --
	Sec Objective Primal, Dual Feasibility Gap Step Path Parameter
	0: -98.8 1 1 1 nan 1
	0: -96.6318 0.432 0.432 0.432 0.578 0.432
	0: -108.421 0.113 0.113 0.113 0.761 0.113
	0: -154.069 0.0292 0.0292 0.0292 0.75 0.0292
	0: -486.968 0.00597 0.00597 0.00597 0.813 0.00597
	0: -991.449 0.0013 0.0013 0.0013 0.803 0.0013
	0: -1223.89 8.92e-05 8.92e-05 8.92e-05 0.963 8.92e-05
	0: -1236.61 1.49e-05 1.49e-05 1.49e-05 0.846 1.49e-05
	0: -1239.65 3.33e-08 3.33e-08 3.33e-08 1 3.33e-08
	Optimization terminated successfully.
	Current function value: -1239.653881
	Iterations: 8
	max memory: 39 M

	--------------------------------------------------------------------------------
	cost: $ 3234 = [610 305 153] * $ [2 4 5.2] lumber, finishing and carpentry
	low medium high [30 40 30] %
	profit $ [1000 1342 1342] average $ 1240 cost $ 3234
	number of desks tables chairs sold:
	[[50 76 76]
	[20 0 0]
	[43 0 0]]
	--------------------------------------------------------------------------------
	cmax: [-2 -4 -5.2 18 12 3 24 16 4 18 12 3]
	x: [610 305 153 50 20 43 76 0.00018 0.00061 76 0.00024 0.00062]
	c*x: [-1220 -1220 -793 900 240 130 1830 0.0028 0.0025 1373 0.0029 0.0019]

	0.48 real 0.86 user 0.28 sys