source: coopr.pysp/trunk/doc/pysp/pyspbody.tex @ 1830

Last change on this file since 1830 was 1830, checked in by wehart, 10 years ago

Update from Dave W.

File size: 57.0 KB
Line 
1% pyspbody.tex
2\def\Argmax{\mathop{\rm argmax}}
3\newcommand{\xBar}{\Vector{\overline{x}}}
4\newcommand{\Bar}[1]{\overline{#1}}
5\newcommand{\Prl}{\mbox{Pr}(\ell)}
6
7\newcommand{\Vector}[1]{\mbox{{\boldmath $#1$}}}
8\newcommand{\Matrix}[1]{\mbox{{\boldmath $#1$}}}
9\newcommand{\VectLen}[1]{\mid\Vector{#1}\mid}
10\newcommand{\Rint}[1]{\int_{#1}^{\infty}}
11%\newcommand{\Size}[1]{\| #1 \|}
12\newcommand{\Size}[1]{\mid #1 \mid}
13\newcommand{\Floor}[1]{\lfloor #1 \rfloor}
14\newcommand{\Ceil}[1]{\left\lceil #1 \right\rceil}
15
16\section{Overview}
17
18The pysp package extends the pyomo modeling language to support multi-stage stochastic programs with
19enumerated scenarios. Pyomo and pysp are Python version 2.6 programs. The underlying algorithm in pysp
20is based on Progressive Hedging (PH) \cite{RockafellarWets}, which progressively computes {\em weights} corresponding to each variable to force convergence. In order to
21specify a program, the user must provide a reference model and a scenario tree.
22
23Provided and the necessary paths have been communicated to the operating system, the command to execute the pysp package is
24of the form:
25\begin{verbatim}
26runph
27\end{verbatim}
28It is possible,
29and generally necessary, to provide command line arguments. The simplest argument causes the program to output help text:
30\begin{verbatim}
31runph --help
32\end{verbatim}
33but notice that there are two dashes before the word ``help.'' Command line arguments are summarized in Section~\ref{cmdargsec}.
34
35\subsection{Reference Model}
36
37The reference model describes the problem for a canonical scenario. It does not make use of, or describe, a scenario index or any information about uncertainty. Typically, it is just the model that would be used if there were only a single scenario. It is given as a pyomo file. Data from an arbitrary scenario is needed to instantiate.
38
39The objective function needs to be separated by stages. The term for each stage should be ``assigned'' (i.e., constrained to be equal to)
40a variable. These variables names are reported in ScenarioStructure.dat so that they can be used for reporting purposes.
41
42\subsection{Scenario Tree}
43
44The scenario tree provides information about the time stages and the nature of the uncertainties. In order to specify a tree, we must
45indicate the time stages at which information becomes available. We also specify the nodes of a tree to indicate which variables
46are associated with which realization at each stage. The data for each scenario is provided in separate data files, one for each scenario.
47
48\section{File Structure}
49
50\begin{itemize}
51\item ReferenceModel.py  (A pyomo model file)
52\item ReferenceModel.dat (data for an arbitrary scenario)
53\item ScenarioStructure.py (do not edit this file)
54\item ScenarioStructure.dat (among other things: the scenario names: Sname)
55\item *Sname.dat   (full data for now) one file for each scenario
56\end{itemize}
57
58In this list we use ``Sname'' as the generic scenario name. The file scenariostructure.dat gives the names of all the scenarios and for each scenario there is a data file with the same name and the suffix ``.dat'' that contains the full specification of data for the scenario.
59
60\subsection{ScenarioStructure.dat}
61
62The file ScenarioStructure.py defines the python sets and parameters needed to describe stochastic elements. This file should not be edited. Data
63to instantiate these sets and parameters is provided by users in the file ScenarioStructure.dat, which can be given in AMPL \cite{ampl} format.
64This file contains the following data:
65
66\begin{itemize}
67\item set Scenarios: List of the names of the scenarios. These names will subsequently be used as indexes in this data file and these names will
68also be used as the root file names for the scenario data files (each of these will have a .dat extension) if the
69parameter ScenarioBasedData is set to True, which is the default.
70\item[]
71item set Stages: List of the names of the time stages, which must be given in time order. In the sequel we will use {\sc StageName} to represent a node name used
72as an index.
73\item[]
74\item set Nodes: List of the names of the nodes in the scenario tree. In the sequel we will use {\sc NodeName} to represent a node name used
75as an index.
76\item[]
77\item param NodeStage: A list of pairs of nodes and stages to indicate the stage for each node.
78\item[]
79\item param Parent: A list of node pairs to indicate the parent of each node that has a parent (the root node will not be listed).
80\item[]
81\item set Children[{\sc NodeName}]: For each node that has children, provide the list of children. No sets will be give for leaf nodes.
82\item[]
83\item param ConditionalProbability: For each node in the scenario tree, give the conditional probability. For the root node it must be given as 1 and for
84the children of any node with children, the conditional probabilities must sum to 1.
85\item[]
86\item param ScenarioLeafNode: A list of scenario and node pairs to indicate the leaf node for each scenario.
87\item[]
88\item set StageVariables[{\sc StageName}]: For each stage, list the pyomo model variables associated with that stage.
89\end{itemize}
90
91The default behavior is one file per scenario and each file has the full data for the scenario. An
92alternative is to specify just the data that changes from the root node in one file per tree node.
93To select this option, add the following line to ScenarioStructure.dat:
94
95\verb|param ScenarioBasedData := False ;|
96
97This will set it up to want a per-node file, something along the lines of what's in \verb|examples/pysp/farmer/NODEDATA|.
98
99\section{Command Line Arguments \label{cmdargsec}}
100
101The basic PH algorithm is controlled by parameters that are set as command line arguments. Note that options begin with a double dash.
102\begin{itemize}
103  \item \verb|-h|, \verb|--help|\\            Show help message and exit.
104  \item \verb|--verbose|\\             Generate verbose output for both initialization and execution.
105  \item \verb|--model-directory|=MODEL\_DIRECTORY\\
106                        The directory in which all model (reference and scenario) definitions are stored. I.e., the ``.py'' files. Default is ".".
107  \item \verb|--instance-directory|=INSTANCE\_DIRECTORY\\
108                        The directory in which all instance (reference and scenario) definitions are stored. I.e., the ``.dat'' files. Default is ".".
109  \item \verb|--solver|=SOLVER\_TYPE\\  The type of solver used to solve scenario sub-problems. Default is cplex.
110  \item \verb|--solver-manager|=SOLVER\_MANAGER\_TYPE\\
111                        The type of solver manager used to coordinate scenario sub-problem solves. Default is serial. This option is changed in parallel applications
112as described in Section~\ref{parallelsec}.
113  \item \verb|--max-iterations|=MAX\_ITERATIONS\\
114                        The maximal number of PH iterations. Default is 100.
115  \item \verb|--default-rho|=DEFAULT\_RHO\\
116                        The default (global) rho for all blended variables. Default is 1.
117  \item \verb|--rho-cfgfile|=RHO\_CFGFILE\\
118                        The name of a configuration script to compute PH rho values. Default is None.
119  \item \verb|--enable-termdiff|-convergence\\
120                        Terminate PH based on the termdiff convergence metric. The convergcne metric is the unscaled sum of differences between variable values and the mean. Default is True.
121  \item \verb|--enable-normalized|-termdiff-convergence\\
122                        Terminate PH based on the normalized termdiff
123                        convergence metric. Each term in the termdiff sum is normalized by the average value (NOTE: it is NOT normalized by the number of scenarios). Default is False.
124  \item \verb|--termdiff-threshold|=TERMDIFF\_THRESHOLD\\
125                        The convergence threshold used in the term-diff and
126                        normalized term-diff convergence criteria. Default is
127                        0.01, which is too low for most problems.
128  \item \verb|--enable-free-discrete-count-convergence|\\
129                        Terminate PH based on the free discrete variable count convergence metric. Default is False.
130  \item \verb|--free-discrete-count-threshold|=FREE\_DISCRETE\_COUNT\_THRESHOLD\\
131                        The convergence threshold used in the criterion based on when the free discrete variable count convergence criterion. Default is 20.
132  \item \verb|--enable-ww-extensions|\\
133                        Enable the Watson-Woodruff PH extensions plugin. Default is False.
134  \item \verb|--ww-extension-cfgfile|=WW\_EXTENSION\_CFGFILE\\
135                        The name of a configuration file for the Watson-Woodruff PH extensions plugin. Default is wwph.cfg.
136  \item \verb|--ww-extension-suffixfile|=WW\_EXTENSION\_SUFFIXFILE\\
137                        The name of a variable suffix file for the Watson-Woodruff PH extensions plugin. Default is wwph.suffixes.
138   \item \verb|--user-defined-extension|=EXTENSIONFILE. Here, "EXTENSIONFILE" is the module name, which is in either the current directory (most likely) or somewhere on your PYTHONPATH. A simple example is "testphextension" plugin that simply prints a message to the screen for each callback. The file testphextension.py can be found in the sources directory and is shown in Section~\ref{ExtensionDetailsSec}. A test of this would be to specify "-user-defined-extension=testphextension", assuming testphextension.py is in your PYTHONPATH or current directory.
139 
140Note that both PH extensions (WW PH and your own) can co-exist; however, the WW plugin will be invoked first.
141 
142  \item \verb|--write-ef|\\            Upon termination, write the extensive form of the model - accounting for all fixed variables.
143  \item \verb|--solve-ef|\\            Following write of the extensive form model, solve it.
144  \item \verb|--ef-output-file|=EF\_OUTPUT\_FILE\\
145                        The name of the extensive form output file (currently only LP format is supported), if writing of the extensive form is enabled. Default is efout.lp.
146  \item \verb|--suppress-continuous-variable-output|\\
147                        Eliminate PH-related output involving continuous variables. Default: no output.
148  \item --keep-solver-files\\   Retain temporary input and output files for scenario sub-problem solves.  Default: files not kept.
149  \item --output-solver-logs\\  Output solver logs during scenario sub-problem solves. Default: no output.
150  \item --output-ef-solver-log\\
151                        Output solver log during the extensive form solve. Default: no output.
152  \item --output-solver-results\\
153                        Output solutions obtained after each scenario sub-problem solve. Default: no output.
154  \item --output-times\\        Output timing statistics for various PH components. Default: no output.
155  \item \verb|--disable-warmstarts|\\
156                  Disable warm-start of scenario sub-problem solves in PH iterations >= 1.
157                  Default=False (i.e., warm starts are the default).
158  \item \verb|--drop-proximal-terms|\\
159                  Eliminate proximal terms (i.e., the quadratic penalty terms) from the weighted PH objective.
160                  Default=False (i.e., but default, the proximal terms are included).
161  \item \verb|--retain-quadratic-binary-terms|\\
162                  Do not linearize PH objective terms involving binary decision variables.
163                  Default=False (i.e., the proximal term for binary variables is linearized by default; this can have some impact on the relaxations during the branch and bound solution process).
164  \item \verb|: --linearize-nonbinary-penalty-terms|\\
165                  Approximate the PH quadratic term for continuous variables with a piece-wise linear function.
166                  Default=False.
167  \item \verb|--checkpoint-interval|\\
168                  The number of iterations between writing of a checkpoint file. Default is 0, indicating never..
169  \item \verb|--restore-from-checkpoint|\\
170                  The name of the checkpoint file from which PH should be initialized. Default is not to restore from a checkpoint.
171  \item \verb|--profile=PROFILE|\\     Enable profiling of Python code.  The value of this option is the number of functions that are summarized. The default is no profiling.
172  \item \verb|--enable-gc|\\           Enable the python garbage collecter. The default is no garbage collection.
173\end{itemize}
174
175\section{Extensions via Callbacks \label{CallbackSec}}
176
177Basic PH can converge slowly, so it is usually advisable to extend it or modify it. In pysp, this is done via the pyomo plug-in
178mechanism. The basic PH implementation provides callbacks that enable access to the data structures used by the algorithm.
179In \S\ref{WWExtensionSec} we describe extensions that are provided with the release. In \S\ref{ExtensionDetailsSec},
180we provide information to power users who may wish to modify or replace the extensions.
181
182\subsection{Watson and Woodruff Extensions \label{WWExtensionSec}}
183
184Watson and Woodruff describe innovations for accelerating PH \cite{phinnovate}, most of which are generalized and
185implemented in the
186file \verb|wwextension.py|, but users generally do not need to know this file name. To invoke the program
187with these additional features, invoke the software with a command of the form:
188\begin{verbatim}
189runph --enable-ww-extensions
190\end{verbatim}
191Many of the examples described in \S\ref{ExampleSec} use this plug-in. The main concept
192is that some integer variables should be fixed as the algorithm progresses for two reasons:
193
194\begin{itemize}
195\item Convergence detection:
196A detailed analysis of PH algorithm behavior on a variety of problem indicates that individual decision variables
197frequently converge to specific,
198fixed values scenarios in early PH iterations. Further, despite interactions among the
199the variables, the value frequently does not change in subsequent PH
200iterations. Such variable ``fixing'' behaviors lead to a potentially powerful, albeit obvious, heuristic: once
201a particular variable has been the same in all scenarios for some number of iterations, fix it to that value.
202For problems where the constraints effectively limit $x$ 
203from both sides, these methods may result in PH encountering infeasible scenario sub-problems even though the
204problem is ultimately feasible.
205\item Cycle detection: When there are integer variables, cycling is sometimes
206encountered, consequently, cycle detection and avoidance mechanisms are required to force
207eventual convergence of the PH algorithm in the mixed-integer case. To detect cycles, we focus on
208repeated occurrences of the weights, implemented using a simple hashing scheme \cite{tabuhash} to
209minimize impact on run-time. Once a cycle in the weight vectors associated with any decision variable is detected,
210the value of that variable is fixed.
211\end{itemize}
212
213Fixing variables aggressively can result in shorter solution times, but can also result in solutions that are not as good. Furthermore, for some problems, aggressive fixing can result in infeasible sub-problems even though the problem is
214ultimately feasible. Many of the parameters discussed in the next subsections control fixing of variables. This is
215discussed in a tutorial in section~\ref{WWTutorialSec}.
216
217\subsubsection{Variable Specific Parameters}
218
219The plug-in makes use of parameters to control behavior at the variable level. Global defaults (to override the defaults stated here) should be set using methods described in \S\ref{ParmSec}. Values for each variable should be set using methods described in \S\ref{SuffixSec}.
220Note that for variable fixing based on convergence detection, iteration zero is treated separately. The parameters
221are as follows:
222
223\begin{itemize}
224\item fix\_continuous\_variables: True or False. If true, fixing applies to all variables. If false, then fixing applies only to discrete variables.
225\item Iter0FixIfConvergedAtLB: 1 (True) or 0 (False). If 1, then discrete variables that are at their lower bound in all scenarios after
226the iteration zero solves will be fixed at that bound.
227\item Iter0FixIfConvergedAtUB: 1 (True) or 0 (False). If 1, then discrete variables that are at their upper bound in all sce<narios after
228the iteration zero solves will be fixed at that bound.
229\item Iter0FixIfConvergedAtNB: = 1  1 (True) or 0 (False). If 1, then discrete variables that are at the same value in all scenarios after
230the iteration zero solves will be fixed at that value, without regard to whether it is a bound. If this is true, it takes precedence. A value of zero, on the other hand, implies that variables will not be fixed at at a non-bound.
231\item FixWhenItersConvergedAtLB:  The number of consecutive PH iterations that discrete variables must be their lower bound in all scenarios before they
232will be fixed at that bound. A value of zero implies that variables will not be fixed at the bound.
233\item FixWhenItersConvergedAtUB: The number of consecutive PH iterations that discrete variables must be their upper bound in all scenarios before they
234will be fixed at that bound. A value of zero implies that variables will not be fixed at the bound.
235\item FixWhenItersConvergedAtNB: The number of consecutive PH iterations that discrete variables must be at the same, consistent value in all scenarios before they
236will be fixed at that value, without regard to whether it is a bound. If this is true, it takes precedence. A value of zero, on the other hand, implies that variables will not be fixed at at a non-bound.
237\item FixWhenItersConvergedContinuous: The number of consecutive PH iterations that continuous variables must be at the same, consistent value in all scenarios before they
238will be fixed at that value. A value of zero implies that continuous variables will not be fixed.
239\item CanSlamToLB: True or False. If True, then slamming can be to the lower bound for any variable.
240\item CanSlamToMin: True or False. If True, then slamming can be to the minimum across scenarios for any variable.
241\item CanSlamToAnywhere: True or False. If True, then slamming can be to any value.
242\item CanSlamToMax: True or False. If True, then slamming can be to the maximum across scenarios for any variable.
243\item CanSlamToUB: True of False. If True, then slamming can be to the upper bound for any variable.
244\end{itemize}
245
246In the event that multiple slam targets are True, then Min and Max trump LB and UB while Anywhere trumps all.
247
248\subsection{General Parameters}
249
250The plug-in also makes use of the following parameters, which should be set using methods described in \S\ref{ParmSec}.
251\begin{itemize}
252\item Iteration0Mipgap: Gives the mipgap to be sent to the solver for iteration zero solves. The default is zero, which causes nothing to be sent to the solver (i.e., the solver uses its default mipgap).
253\item InitalMipGap: Gives the mipgap to be sent to the solver for iteration one solves. The default is zero, which causes nothing to be sent to the solver (i.e., the solver uses its default mipgap). If not zero, then this gap will be used as the starting point for the mipgap to change on each PH iteration in proportion to the convergence termination criterion so that when the criterion for termination is met
254the mipgap will be at (or near) the parameter value of FinalMipGap. If the InitialMipGap is significantly higher than the Iteration0MipGap parameter, the PH algorithm may
255perform poorly. This is because the iteration k-1 solutions are used to warm start iteration k-1 solves. If the iteration 1 mipgap is much higher than the iteration 0 mipgap, the iteration zero
256solution, although not optimal for the iteration one objective, might be within the mipgap. Default: 0.10.
257\item FinalMipGap: The target for the mipgap when PH has converged. Default: 0.001.
258\item PH\_Iters\_Between\_Cycle\_Slams: controls the number of iterations to wait after a variable is slammed due to
259hash hits that suggest convergence. Zero indicates unlimited slams per cycle. Default: 1.
260\item SlamAfterIter: Iteration number after which one variable every other iteration will be slammed to force convergence. Default: the number of scenarios.
261\item hash\_hit\_len\_to\_slam: Ignore possible cycles for which the only evidence of a cycle is less than this. Also, ignore cycles if any variables have been fixed in the previous hash\_hit\_len\_to\_slam PH iterations. Default:
262the number of scenarios. This default is often not a good choice. For many problems with a lot of scenarios,
263a value like 10 or 20 might be a lot better if rapid convergence is desired.
264\item W\_hash\_history\_len: This obscure parameter controls how far back the code will look to see if there is a possible cycle. Default: max(100, number of scenarios).
265\end{itemize}
266
267\subsubsection{Setting Parameter Values \label{ParmSec}}
268
269The parameters of PH and of any callbacks can be changed using the file wwph.cfg, which is executed by the python interpreter after PH has initialized. This is
270a potentially powerful and/or dangerous procedure because it gives an opportunity
271to change anything using the full features of python and pyomo.
272
273\subsubsection{Setting Suffix Values \label{SuffixSec}}
274
275Suffixes are set using the data file named \verb|wwph.suffixes| using this syntax:
276
277\begin{verbatim}
278VARSPEC SUFFIX VALUE
279\end{verbatim}
280
281where VARSPEC is replaced by a variable specification, SUFFIX is replaced by a suffix name and VALUE is replaced by the value of the suffix for
282that variable or those variables. Here is an example:
283
284\begin{verbatim}
285   Delta CanSlamToLB False
286   Gamma[*,Ano1] SlammingPriority 10
287   Gamma[*,Ano2] SlammingPriority 20
288   ...
289\end{verbatim}
290
291
292\subsection{Callback Details \label{ExtensionDetailsSec}}
293
294Most users of pysp can skip this subsection.
295A callback class definition named iphextension is in the file iphextension.py and can be used to implement callbacks at a
296variety of points in PH. For example, the method post\_iteration\_0\_solves is called immediately after all iteration zero
297solves, but before averages and weights have been computed while the method post\_iteration\_0 is called after averages and
298weights based on iteration zero have been computed. The file iphextension is in the coopr/pysp directory and is not
299intended to be edited by users.
300
301The user defines a class derived from SingletonPlugin that implements iphextension. Its name is given to ph as an option. This class will be automatically instantiated by ph.
302It has access to data and methods in the PH class, which are defined in the file ph.py. An example of such a class is in the
303file named testphextension.py in the pysp example directory.
304
305A user defined extension file can be incorporated by using the command line option:
306\verb|--user-defined-extension=EXTENSIONFILE|. Here, "EXTENSIONFILE" is the module name, which is in either the current directory (most likely) or somewhere on your PYTHONPATH. A simple example is "testphextension" plugin that simply prints a message to the screen for each callback. The file testphextension.py can be found in the sources directory and given in Section~\ref{callbacksec}. An easy test of this would be to specify "-user-defined-extension=testphextension" and you should
307note the the ``.py'' file extension is not included on the runph command line.
308 
309Both your own extension and the WWPH extensions can co-exist; however, the WW plugin will be invoked first at each callback point if both are included.
310
311Here are the callbacks:
312\begin{itemize}
313\item post\_iteration\_0\_solves: Called after iteration zero solutions and some statistics such as averages have been computed, but before weights are updated.
314\item post\_iteration\_0: Called after all processing for iteration zero is complete.
315\item post\_iteration\_k\_solves: Called after solutions some statistics such as averages have been computed, but before weights are updated for iterations after iteration zero.
316\item post\_iteration\_k: Called after all processing for each iteration after iteration 0 is complete.
317\item post\_ph\_execution: Called execution is complete.
318\end{itemize}
319
320Users interested in writing their own extensions will probably want to refer to the source file ph.py to get a deeper understanding of when the callback occur.
321
322\section{Examples \label{ExampleSec}}
323
324A number of examples are provided with pysp.
325
326\subsection{Farmer Example}
327
328This two-stage example is composed of models and data for the "Farmer" stochastic program, introduced in Section 1.1 of
329"Introduction to Stochastic Programming" by Birge and
330Louveaux \cite{spbook}.
331
332\begin{itemize}
333\item ReferenceModel.py: a single-scenario model for the SP
334\item ReferenceModel.dat: a single-scenario data file for the SP (any scenario will do - used to flush out variable and constraint index sets)
335\item ScenarioStructure.py: defines the scenario tree structure for the SP. SHOULD NOT BE MODIFIED.
336\item ScenarioStructure.dat: data file defining the scenario tree.
337\item AboveAverageScenario.dat: one of the scenario data files.
338\item BelowAverageScenario.dat: one of the scenario data files.
339\item AverageScenario.dat: one of the scenario data files.
340\end{itemize}
341
342The command runph executes PH, assuming the ReferenceModel.* and ScenarioStructure.* files are present and correct.
343This example is probably in a directory with a name something like:
344
345\begin{verbatim}
346pyomo\examples\pysp\farmer
347\end{verbatim}
348where the directory \verb|pyomo| might be \verb|pyomodist| or something else depending on how you installed pyomo and pysp.
349
350The data is in a subdirectory called SCENARIODATA.
351
352To
353invoke PH for this problem, connect to this farmer directory and use the command:
354
355\begin{verbatim}
356runph --model-directory=models --instance-directory=scenariodata
357\end{verbatim}
358
359\subsection{Sizes Example}
360
361This two-stage example is composed of models and data for the "Sizes" stochastic program \cite{sizes,lokwood}, which consists of the following files:
362
363\begin{itemize}
364\item phdriver.ph: executes ph. Do not modify this file unless you an algorithm developer and you
365need to edit a line to speficy a plug-in that you have written. It is presently configured to use the WWPHExtensions.py plug-in the implements the Watson and Woodruff extensions.
366\item wwph.cfg: replace default algorithm parameter values for the Watson and Woodruff extensions.
367\item wwph.suffixes: sets algorithm parameter values at the variables level for the Watson and Woodruff extensions.
368\item ReferenceModel.py: a single-scenario model for the SP
369\item ReferenceModel.dat: a single-scenario data file for the SP (any scenario will do - used to flush out variable and constraint index sets)
370\item ScenarioStructure.py: defines the scenario tree structure for the SP. SHOULD NOT BE MODIFIED.
371\item ScenarioStructure.dat: data file defining the scenario tree.
372\item Scenario1.dat: one of the scenario data files.
373\item Scenario2.dat: one of the scenario data files.
374\item ...
375\end{itemize}
376
377The file phdriver.py executes PH, assuming the ReferenceModel.* and ScenarioStructure.* files are present and correct.
378This example is probably in a directory with a name something like:
379
380\begin{verbatim}
381pyomo\packages\coopr\examples\pysp\sizes
382\end{verbatim}
383where the directory \verb|pyomo| might be \verb|pyomodist| or something else depending on how you installed pyomo and pysp.
384
385The data for a three scenario version is in a subdirectory called SIZES3 and a ten scenario
386dataset is in SIZES10.
387
388To invoke PH for the 10 scenario problem, connect to the sizes directory and use the command:
389
390\begin{verbatim}
391runph --model-directory=models --instance-directory=SIZES10
392\end{verbatim}
393
394
395\subsection{Forestry Example}
396
397This four-stage example is composed of models and data for the ``forestry'' stochastic program \cite{}, which consists of the following files:
398
399\begin{itemize}
400\item phdriver.ph: executes ph. Do not modify this file unless you an algorithm developer and you
401need to edit a line to specify a plug-in that you have written. It is presently configured to use the WWPHExtensions.py plug-in the implements the Watson and Woodruff extensions.
402\item wwph.cfg: replace default algorithm parameter values for the Watson and Woodruff extensions.
403\item wwph.suffixes: sets algorithm parameter values at the variables level for the Watson and Woodruff extensions.
404\item ReferenceModel.py: a single-scenario model for the SP
405\item ReferenceModel.dat: a single-scenario data file for the SP (any scenario will do - used to flush out variable and constraint index sets)
406\item ScenarioStructure.py: defines the scenario tree structure for the SP. SHOULD NOT BE MODIFIED.
407\item ScenarioStructure.dat: data file defining the scenario tree.
408\item Scenario1.dat: one of the scenario data files.
409\item Scenario2.dat: one of the scenario data files.
410\item ...
411\end{itemize}
412
413The file phdriver.py executes PH, assuming the ReferenceModel.* and ScenarioStructure.* files are present and correct.
414This example is probably in a directory with a name something like:
415
416\begin{verbatim}
417pyomo\packages\coopr\examples\pysp\forestry
418\end{verbatim}
419where the directory \verb|pyomo| might be \verb|pyomodist| or something else depending on how you installed pyomo and pysp.
420
421There are two families of instances: ``Chile'' and ``Davis,'' each with four stages and eighteen scenarios. This is also a small
422two-stage, four scenario instances in the subdirectory DAVIS2STAGE.
423
424This full 18 scenario problem instance takes too long without the Watson Woodruff extensions, but to
425invoke PH for this problem without them, connect to the forestry directory and use the command:
426
427\begin{verbatim}
428 runph --models-directory=models --instancs-directory=chile
429\end{verbatim}
430and to run with the extensions, use
431\begin{verbatim}
432runph --model-directory=models --instance-directory=davis \
433--enable-ww-extensions --ww-extension-cfgfile=config/wwph.cfg \
434--ww-extension-suffixfile=config/wwph.suffixes
435\end{verbatim}
436
437
438\section{Tutorial: Parameters for Watson and Woodruff PH Extensions \label{WWTutorialSec}}
439
440The parameters for the PH extensions in WWPHExtensions.py provide the user with considerable control over
441how and under what conditions variables are fixed during the PH algorithm execution. Often, some variables converge
442to consistent values during early iterations and can be fixed at these values without affecting quality very much and without
443affecting feasibility at all. Also, the algorithm may need to fix some variables during execution in order to break cycles.
444In both cases, guidance from the user concerning which classes of variables can and/or should be fixed under various circumstances can be very helpful.
445
446The overarching goal is to support industrial and government users who solve the same model repeatedly with different, but
447somewhat similar, data each time. In such settings, it behooves the modeler to consider tradeoffs between speed, solution quality
448and feasibility and create at least one good set of directives and parameters for the PH algorithm. In some cases, a user may want more than
449one set of directives and parameters depending on whether speed of execution or quality of solution are more important. Iteration zero is controlled separately because often the absence of the quadratic term results in faster solves for this
450iteration and fixing variables after the iteration has the maximum possible impact on speedup.
451
452In order to discuss these issues, we consider an example.
453
454\subsection{Sizes Example}
455
456The Sizes example is simple and small. In fact, the instances that we will consider are so small that there is really
457no need to use the PH algorithm since the extensive form can be solved directly in a few minutes. However, it serves as a
458vehicle for introducing the concepts.
459
460This description and formulation comes from the original paper by Jorjani, Scott and Woodruff \cite{sizes}.
461
462If demand constraints for the current period are
463based on firm orders but future demands are based on forecasts or
464conjecture, then they should not be treated in the same fashion. We
465must recognize that future demand constraints are stochastic. That is
466to say that they should be modeled as random variables.
467It may not be reasonable or useful to consider the entire demand probability
468distribution functions. It may not be reasonable because there may not be
469sufficient data to estimate an entire distribution. It may not be useful
470because the essence of the stochastics may be captured by specifying a
471small number of representative {\em scenarios}. We assume that scenarios
472are specified by giving a full set of random variable realizations and a
473corresponding probability.
474We index the scenario set, ${\cal L}$, by $\ell$ and
475refer to the probability of occurrence of $\ell$ (or, more accurately, a
476realization ``near'' scenario $\ell$) as $\Prl$.
477We refer to solution systems that satisfy constraints with probability one as
478{\em admissible}.
479
480In addition to modeling stochastics, we would like to model {\em
481recourse} as well. That is, we would like to model the ability of
482decision makers to make use of new information (e.g., orders) at the
483start of each planning period.  We allow our solution vectors to
484depend on the scenario that is realized, but we do not want to assume
485prescience. We refer to a system of solution vectors as {\em
486implementable} if for all decision times $t$, the solution vector
487elements corresponding to period $1,\ldots,t$ are constant with
488respect to information that becomes available only after stage $t$.
489We refer to the set of implementable solutions as ${\cal N}_{\cal L}$.  It is
490possible to require implementable solutions by adding {\em
491non-anticipatitivity constraints}, but we will instead make use of
492solution procedures that implicitly guarantee implementable solutions.
493
494In the stochastic, multi-period formulation that follows the objective
495is to minimize expected costs.  We invoke the network equivalence
496given earlier to drop the explicit requirement that $\Vector{x}$ and
497$\Vector{y}$ be integers. Variables and data are subscripted with a
498period index $t$ that takes on values up to $T$. To model the idea
499that sleeves produced in one period can be used as-is or cut in
500subsequent periods, we use $x_{ijt}$ to indicate that sleeves of
501length index $i$ are to be used without cutting in period $t$ if $i=j$
502and with cutting otherwise.  The $\Vector{y}$ vector gives production
503quantities for each length in each period without regard to the period
504in which they will be used (and perhaps cut). The formulation is
505essentially an extension of ILP except that a capacity constraint
506must be added in the multiple period formulation. Holding costs could
507be added, but an additional subscript becomes necessary without the
508benefit of any additional insight. As an aside, note that the addition
509of holding costs would add a large number of continuous variables, but
510no new integers so the impact on computational performance would not
511be catastrophic.
512
513\[
514\min\sum_{\ell\in{\cal L}}\Prl
515     \sum_{t}^{T}\left[\sum_{i=1}^{N}\left(sz_{it\ell}+p_{i}y_{it\ell}\right)
516    +r\sum_{j<i}x_{ijt\ell}\right] \;\;\;\;\mbox{(SMIP)}
517\]
518subject to
519\begin{eqnarray}
520\sum_{j \geq i}x_{ijt\ell} & \geq & d_{it\ell} 
521          \;\;\;\; \ell \in {\cal L}, \; i=1,\ldots,N, \; t=1,\ldots,T \label{Dconstr} \\
522\sum_{t'\leq t}\left[\sum_{j \leq i}x_{ijt'\ell} - y_{it'\ell}\right] & \leq & 0
523          \;\;\;\; \ell \in {\cal L}, \; i=1,\ldots,N, \; t=1,\ldots,T \label{Pconstr} \\
524y_{it\ell} - Mz_{it\ell} & \leq & 0
525          \;\;\;\; \ell \in {\cal L}, \;i=1,\ldots,N, \; t=1,\ldots,T \label{Mconstr} \\
526\sum_{i=1}^{N}y_{it\ell} & \leq & c_{t\ell}
527          \;\;\;\; \ell \in {\cal L}, \; t=1,\ldots,T \label{Cconstr} \\
528z_{it\ell} & \in & \{0,1\} 
529          \;\;\;\; \ell \in {\cal L}, \; i=1,\ldots,N, \; t=1,\ldots,T \label{intconstr} \\
530\Vector{x}, \Vector{y}, \Vector{z} & \in & {\cal N}_{\cal L} 
531\end{eqnarray}
532
533Bear in mind that solution vector elements corresponding
534to periods two through $T$ are not actually intended for use, they
535are computed just to see the effect that period 1 (the current period)
536decision would have on future optimal behavior. At the start of period 2
537-- or at any other time -- the decision maker would run the model again
538with updated demands for the current period and new scenario estimates.
539
540\subsection{ReferenceModel.py}
541
542Here is the single scenario reference model:
543
544{\small
545\begin{verbatim}
546#
547# This is the two-period version of the SIZES optimization model.
548#
549
550from coopr.pyomo import *
551
552#
553# Model
554#
555
556model = Model()
557
558#
559# Parameters
560#
561
562# the number of product sizes.
563model.NumSizes = Param(within=NonNegativeIntegers)
564
565def product_sizes_rule(model):
566    ans = set()
567    for i in range(1, model.NumSizes()+1):
568        ans.add(i)
569    return ans
570
571# the set of sizes, labeled 1 through NumSizes.
572model.ProductSizes = Set(initialize=product_sizes_rule)
573
574# the deterministic demands for product at each size.
575model.DemandsFirstStage = Param(model.ProductSizes, within=NonNegativeIntegers)
576model.DemandsSecondStage = Param(model.ProductSizes, within=NonNegativeIntegers)
577
578# the unit production cost at each size.
579model.UnitProductionCosts = Param(model.ProductSizes, within=NonNegativeReals)
580
581# the setup cost for producing any units of size i.
582model.SetupCosts = Param(model.ProductSizes, within=NonNegativeReals)
583
584# the unit penalty cost of meeting demand for size j with larger size i.
585model.UnitPenaltyCosts = Param(model.ProductSizes, within=NonNegativeReals)
586
587# the cost to reduce a unit i to a lower unit j.
588model.UnitReductionCost = Param(within=NonNegativeReals)
589
590# a cap on the overall production within any time stage.
591model.Capacity = Param(within=PositiveReals)
592
593# a derived set to constrain the NumUnitsCut variable domain.
594def num_units_cut_domain_rule(model):
595   ans = set()
596   for i in range(1,model.NumSizes()+1):
597      for j in range(1, i+1):
598         ans.add((i,j))   
599   return ans
600
601model.NumUnitsCutDomain = Set(initialize=num_units_cut_domain_rule, dimen=2)
602
603#
604# Variables
605#
606
607# are any products at size i produced?
608model.ProduceSizeFirstStage = Var(model.ProductSizes, domain=Boolean)
609model.ProduceSizeSecondStage = Var(model.ProductSizes, domain=Boolean)
610
611# NOTE: The following (num-produced and num-cut) variables are implicitly integer
612#       under the normal cost objective, but with the PH cost objective, this isn't
613#       the case.
614
615# the number of units at each size produced.
616model.NumProducedFirstStage = Var(model.ProductSizes, domain=NonNegativeIntegers)
617model.NumProducedSecondStage = Var(model.ProductSizes, domain=NonNegativeIntegers)
618
619# the number of units of size i cut (down) to meet demand for units of size j.
620model.NumUnitsCutFirstStage = Var(model.NumUnitsCutDomain, domain=NonNegativeIntegers)
621model.NumUnitsCutSecondStage = Var(model.NumUnitsCutDomain, domain=NonNegativeIntegers)
622
623# stage-specific cost variables for use in the pysp scenario tree / analysis.
624model.FirstStageCost = Var(domain=NonNegativeReals)
625model.SecondStageCost = Var(domain=NonNegativeReals)
626
627#
628# Constraints
629#
630
631# ensure that demand is satisfied in each time stage, accounting for cut-downs.
632def demand_satisfied_first_stage_rule(i, model):
633   return (0.0, \
634           sum(model.NumUnitsCutFirstStage[j,i] \
635               for j in model.ProductSizes if j >= i) - model.DemandsFirstStage[i], \
636           None)
637
638def demand_satisfied_second_stage_rule(i, model):
639   return (0.0, \
640           sum(model.NumUnitsCutSecondStage[j,i] \
641               for j in model.ProductSizes if j >= i) - model.DemandsSecondStage[i], \
642           None)
643
644model.DemandSatisfiedFirstStage = \
645                Constraint(model.ProductSizes, rule=demand_satisfied_first_stage_rule)
646model.DemandSatisfiedSecondStage = \
647                Constraint(model.ProductSizes, rule=demand_satisfied_second_stage_rule)
648
649# ensure that you don't produce any units if the decision has been made to disable production.
650def enforce_production_first_stage_rule(i, model):
651   # The production capacity per time stage serves as a simple upper bound for "M".
652   return (None, \
653           model.NumProducedFirstStage[i] - model.Capacity * model.ProduceSizeFirstStage[i], \
654           0.0)
655
656def enforce_production_second_stage_rule(i, model):
657   # The production capacity per time stage serves as a simple upper bound for "M".   
658   return (None, \
659           model.NumProducedSecondStage[i] - model.Capacity * model.ProduceSizeSecondStage[i], \
660           0.0)
661
662model.EnforceProductionBinaryFirstStage = \
663                  Constraint(model.ProductSizes, rule=enforce_production_first_stage_rule)
664model.EnforceProductionBinarySecondStage = \
665                  Constraint(model.ProductSizes, rule=enforce_production_second_stage_rule)
666
667# ensure that the production capacity is not exceeded for each time stage.
668def enforce_capacity_first_stage_rule(model):
669   return (None, \
670           sum(model.NumProducedFirstStage[i] for i in model.ProductSizes) - model.Capacity, \
671           0.0)
672
673def enforce_capacity_second_stage_rule(model):
674   return (None, \
675           sum(model.NumProducedSecondStage[i] for i in model.ProductSizes) - model.Capacity, \
676           0.0)   
677
678model.EnforceCapacityLimitFirstStage = Constraint(rule=enforce_capacity_first_stage_rule)
679model.EnforceCapacityLimitSecondStage = Constraint(rule=enforce_capacity_second_stage_rule)
680
681# ensure that you can't generate inventory out of thin air.
682def enforce_inventory_first_stage_rule(i, model):
683   return (None, \
684           sum(model.NumUnitsCutFirstStage[i,j] \
685               for j in model.ProductSizes if j <= i) - model.NumProducedFirstStage[i], \
686           0.0)
687
688def enforce_inventory_second_stage_rule(i, model):
689   return (None, \
690           sum(model.NumUnitsCutFirstStage[i,j] \
691               for j in model.ProductSizes \
692                  if j <= i) + sum(model.NumUnitsCutSecondStage[i,j] \
693                     for j in model.ProductSizes if j <= i) \
694               - model.NumProducedFirstStage[i] - model.NumProducedSecondStage[i], \
695           0.0)
696
697model.EnforceInventoryFirstStage = \
698               Constraint(model.ProductSizes, rule=enforce_inventory_first_stage_rule)
699model.EnforceInventorySecondStage = \
700               Constraint(model.ProductSizes, rule=enforce_inventory_second_stage_rule)
701
702# stage-specific cost computations.
703def first_stage_cost_rule(model):
704   production_costs = \
705            sum(model.SetupCosts[i] * model.ProduceSizeFirstStage[i] \
706                + model.UnitProductionCosts[i] * model.NumProducedFirstStage[i] \
707                for i in model.ProductSizes)
708   cut_costs = \
709             sum(model.UnitReductionCost * model.NumUnitsCutFirstStage[i,j] \
710                 for (i,j) in model.NumUnitsCutDomain if i != j)
711   return (model.FirstStageCost - production_costs - cut_costs) == 0.0
712
713model.ComputeFirstStageCost = Constraint(rule=first_stage_cost_rule)
714
715def second_stage_cost_rule(model):
716   production_costs = \
717                sum(model.SetupCosts[i] * model.ProduceSizeSecondStage[i] \
718                    + model.UnitProductionCosts[i] * model.NumProducedSecondStage[i] \
719                    for i in model.ProductSizes)
720   cut_costs = \
721             sum(model.UnitReductionCost * model.NumUnitsCutSecondStage[i,j] \
722                 for (i,j) in model.NumUnitsCutDomain if i != j)
723   return (model.SecondStageCost - production_costs - cut_costs) == 0.0   
724
725model.ComputeSecondStageCost = Constraint(rule=second_stage_cost_rule)
726
727#
728# Objective
729#
730
731def total_cost_rule(model):
732   return (model.FirstStageCost + model.SecondStageCost)
733
734model.TotalCostObjective = Objective(rule = total_cost_rule, sense=minimize)
735\end{verbatim}
736}
737
738\subsection{ReferenceModel.dat}
739
740This file establishes the data for a representative instance. The main thing to notice is that the indexes for
741sizes happen to be given as integers, which is not required: they could have been strings.
742
743{\small
744\begin{verbatim}
745#ReferenceModel.dat
746param NumSizes := 10 ;
747
748param Capacity := 200000 ;
749
750param DemandsFirstStage := 1 2500 2 7500 3 12500 4 10000 5 35000
751                           6 25000 7 15000 8 12500 9 12500 10 5000 ;
752
753param DemandsSecondStage := 1 2500 2 7500 3 12500 4 10000 5 35000
754                            6 25000 7 15000 8 12500 9 12500 10 5000 ;
755
756param UnitProductionCosts := 1 0.748 2 0.7584 3 0.7688 4 0.7792 5 0.7896
757                             6 0.8 7 0.8104 8 0.8208 9 0.8312 10 0.8416 ;
758
759param SetupCosts := 1 453 2 453 3 453 4 453 5 453 6 453 7 453 8 453 9 453 10 453 ;
760
761param UnitReductionCost := 0.008 ;
762\end{verbatim}
763
764\subsection{ScenarioStucture.dat}
765
766Here is the data file that describes the stochastics:
767
768\begin{verbatim}
769# IMPORTANT - THE STAGES ARE ASSUMED TO BE IN TIME-ORDER.
770
771set Stages := FirstStage SecondStage ;
772
773set Nodes := RootNode
774             Scenario1Node
775             Scenario2Node
776             Scenario3Node
777             Scenario4Node
778             Scenario5Node
779             Scenario6Node
780             Scenario7Node
781             Scenario8Node
782             Scenario9Node
783             Scenario10Node ;
784
785param NodeStage := RootNode             FirstStage
786                   Scenario1Node        SecondStage
787                   Scenario2Node        SecondStage
788                   Scenario3Node        SecondStage
789                   Scenario4Node        SecondStage
790                   Scenario5Node        SecondStage
791                   Scenario6Node        SecondStage
792                   Scenario7Node        SecondStage
793                   Scenario8Node        SecondStage
794                   Scenario9Node        SecondStage
795                   Scenario10Node       SecondStage ;
796
797set Children[RootNode] := Scenario1Node
798                          Scenario2Node
799                          Scenario3Node
800                          Scenario4Node
801                          Scenario5Node
802                          Scenario6Node
803                          Scenario7Node
804                          Scenario8Node
805                          Scenario9Node
806                          Scenario10Node ;
807
808param ConditionalProbability := RootNode          1.0
809                                Scenario1Node     0.10
810                                Scenario2Node     0.10
811                                Scenario3Node     0.10
812                                Scenario4Node     0.10
813                                Scenario5Node     0.10
814                                Scenario6Node     0.10
815                                Scenario7Node     0.10
816                                Scenario8Node     0.10
817                                Scenario9Node     0.10
818                                Scenario10Node    0.10 ;
819
820set Scenarios := Scenario1
821                 Scenario2
822                 Scenario3
823                 Scenario4
824                 Scenario5
825                 Scenario6
826                 Scenario7
827                 Scenario8
828                 Scenario9
829                 Scenario10 ;
830
831param ScenarioLeafNode := Scenario1  Scenario1Node 
832                          Scenario2  Scenario2Node 
833                          Scenario3  Scenario3Node 
834                          Scenario4  Scenario4Node 
835                          Scenario5  Scenario5Node 
836                          Scenario6  Scenario6Node 
837                          Scenario7  Scenario7Node 
838                          Scenario8  Scenario8Node 
839                          Scenario9  Scenario9Node 
840                          Scenario10 Scenario10Node ;
841
842set StageVariables[FirstStage] := ProduceSizeFirstStage[*]
843                                  NumProducedFirstStage[*]
844                                  NumUnitsCutFirstStage[*,*] ;
845set StageVariables[SecondStage] := ProduceSizeSecondStage[*]
846                                   NumProducedSecondStage[*]
847                                   NumUnitsCutSecondStage[*,*] ;
848
849param StageCostVariable := FirstStage FirstStageCost
850                           SecondStage SecondStageCost ;
851                                   
852\end{verbatim}
853}
854
855\subsection{wwph.cfg}
856
857This file provides default values. Even though most of them will be overridden at the variable level, it makes sense
858to provide instance specific defaults. For this problem, it does not seem helpful or prudent to do any fixing of continuous
859variables since they only used to compute the objective function terms. The \verb|ProduceSizeFirstStage| variables are binary
860and fixing their values would seem to determine the value of the other values. The other first stage variables are general integers. The \verb|ProduceSizeFirstStage| can safely be fixed because provided that the largest size is not fixed at zero, there is little risk of infeasibility since larger sizes can be cut (at a cost) to meet demand for smaller sizes. Consequently, it is safe to let the algorithm slam those
861variables as needed. Slamming the other variables is riskier because they could get slammed to values that don't make
862sense given the values of the \verb|ProduceSizeFirstStage| variables.
863
864Fixing variables that have converged will speed the algorithm, perhaps resulting in a less desirable
865objective function value. It would make sense to tend to fix the \verb|ProduceSizeFirstStage| before fixing the others to avoid
866inconsistencies and because the The \verb|ProduceSizeFirstStage| variables have a strong tendency to imply values for the other
867variables.
868
869Here is a sensible and internally consistent, if a bit conservative, wwph.cfg file:
870
871{\small
872\begin{verbatim}
873# wwph.cfg
874# python  commands that are executed by the wwphextensions.py file
875# to set parameters and parameter defaults
876
877self.fix_continuous_variables = False
878
879self.Iteration0MipGap = 0.02
880self.InitialMipGap = 0.10
881self.FinalMipGap = 0.001
882
883# for all six of these, zero means don't do it.
884self.Iter0FixIfConvergedAtLB = 0 # 1 or 0
885self.Iter0FixIfConvergedAtUB = 0  # 1 or 0
886self.Iter0FixIfConvergedAtNB = 0  # 1 or 0 (converged to a non-bound)
887self.FixWhenItersConvergedAtLB = 0
888self.FixWhenItersConvergedAtUB = 25
889self.FixWhenItersConvergedAtNB = 0  # converged to a non-bound
890self.FixWhenItersConvergedContinuous = 0
891     
892# "default" slamming parms
893# True and False are the options (case sensitive)
894self.CanSlamToLB = False
895self.CanSlamToMin = False
896self.CanSlamToAnywhere = False
897self.CanSlamToMax = True
898self.CanSlamToUB = False
899self.PH_Iters_Between_Cycle_Slams = 5
900
901# the next line will try to force at least one variable to be
902# fixed every other iteration after iteration 50
903# if anything can be slammed
904self.SlamAfterIter = 50
905
906self.hash_hit_len_to_slam = 50
907self.W_hash_history_len = 100
908\end{verbatim}
909}
910
911There are a number of things to notice about the contents of this file. Since it is executed
912as python code, the syntax matters. Users should changes values of numbers of change True to False, but
913everything else should not be edited with the exception of comments, which is any text after a sharp sign
914(sometimes called a pound sign). Changes to this file should be tested incrementally because errors are trapped by
915the python interpreter and may be difficult for non-python programmers to decipher.
916
917This particular example file allows variables to be fixed if all scenarios have agreed on the
918upper bound for five iterations in a row. Since the \verb|ProduceSizeFirstStage| variables are the only discrete
919variables in the first stage with an upper bound, they are the only variables affected by this.
920
921This example allows slamming, but only to the max across scenarios. This is different than the upper bound, even for
922binary variables, because a variable could be selected for slamming for which all scenarios have agreed on the same value
923(which could be the lower lower bound). Data providing an override for this default as well as control over the selection priority for variables to slam are provided in the wwph.suffixes file.
924
925\subsection{wwph.suffixes}
926
927{\small
928\begin{verbatim}
929# Optional suffixes to help control PH
930
931# Text triples specifying (variable/variable-index, suffix-name, suffix-value) tuples.
932
933# override the defaults given in wwph.cfg as needed
934# If no scenario needs the smallest size to be produced, then just forget it
935ProduceSizeFirstStage[1] Iter0FixIfConvergedAtLB 1
936
937# if all scenarios need either of the two largest sizes, just lock them in
938ProduceSizeFirstStage[9] Iter0FixIfConvergedAtUB 1
939ProduceSizeFirstStage[10] Iter0FixIfConvergedAtUB 1
940
941# and be aggressive with the smallest size after iteration 0
942ProduceSizeFirstStage[1] FixWhenItersConvergedAtLB 8
943
944# if the production quantities have been fixed a long time, fix them
945NumProducedFirstStage[*] FixWhenItersConvergedAtNB 20
946
947# don't allow slamming of variables other than ProduceSize
948# for completeness, disallow all slamming destinations
949
950NumProducedFirstStage[*] CanSlamToLB False
951NumProducedFirstStage[*] CanSlamToMin False
952NumProducedFirstStage[*] CanSlamToAnywhere False
953NumProducedFirstStage[*] CanSlamToMax False
954NumProducedFirstStage[*] CanSlamToUB False
955
956NumUnitsCutFirstStage[*,*] CanSlamToLB False
957NumUnitsCutFirstStage[*,*] CanSlamToMin False
958NumUnitsCutFirstStage[*,*] CanSlamToAnywhere False
959NumUnitsCutFirstStage[*,*] CanSlamToMax False
960NumUnitsCutFirstStage[*,*] CanSlamToUB False
961
962# higher priority means slam these first
963ProduceSizeFirstStage[1] SlammingPriority 1
964ProduceSizeFirstStage[2] SlammingPriority 2
965ProduceSizeFirstStage[3] SlammingPriority 3
966ProduceSizeFirstStage[4] SlammingPriority 4
967ProduceSizeFirstStage[5] SlammingPriority 5
968ProduceSizeFirstStage[6] SlammingPriority 6
969ProduceSizeFirstStage[7] SlammingPriority 7
970ProduceSizeFirstStage[8] SlammingPriority 8
971ProduceSizeFirstStage[9] SlammingPriority 9
972ProduceSizeFirstStage[10] SlammingPriority 10
973
974\end{verbatim}
975}
976
977The first thing to notice is that this is a data file and not a file of Python comments.
978Consequently, simpler messages are provided if there are errors, but the file can be verbose
979and a computer program or spreadsheet might be required to produce this data file. The next thing to notice
980is that indexes can be specified either using valid index values, or wildcards.
981
982This file overrides the defaults to allow some fixing after iteration zero. Fixing the smallest
983size (index 1) to zero cannot result in infeasibility because larger sizes can be cut to
984satisfy demand for that size. Along the same lines, aggressively fixing the
985production indicator for larger sizes to 1 is also
986safe and perhaps not sub-optimal if all scenarios ``want'' those sizes anyway.
987
988\subsection{rhosetter.cfg}
989
990\begin{verbatim}
991# rhosetter.cfg
992
993# users probably want this line intact so they can use model_instance
994model_instance = self._model_instance
995
996MyRhoFactor = 0.1
997
998for i in model_instance.ProductSizes:
999   self.setRhoAllScenarios(model_instance.ProduceSizeFirstStage[i], model_instance.SetupCosts[i] * MyRhoFactor)
1000   self.setRhoAllScenarios(model_instance.NumProducedFirstStage[i], model_instance.UnitProductionCosts[i] * MyRhoFactor)
1001   for j in model_instance.ProductSizes:
1002      if j <= i:
1003         self.setRhoAllScenarios(model_instance.NumUnitsCutFirstStage[i,j], model_instance.UnitReductionCost * MyRhoFactor)
1004
1005
1006\end{verbatim}
1007
1008\subsection{Execution}
1009
1010To run this example, connect to the sizes directory, which is something like:
1011\begin{verbatim}
1012pyomo\packages\coopr\examples\pysp\sizes
1013\end{verbatim}
1014where the directory \verb|pyomo| might be \verb|pyomodist| or something else depending on how you installed pyomo and pysp.
1015
1016Then use
1017\begin{verbatim}
1018runph --model-directory=models --instance-directory=SIZES10 \
1019--enable-ww-extensions --ww-extension-cfgfile=config/wwph.cfg \
1020--ww-extension-suffixfile=config/wwph.suffixes \
1021--rho-cfgfile=config/rhosetter.cfg
1022\end{verbatim}
1023
1024Since this instance is so small by modern standards, the enhancement are not needed and may even increase the total solution
1025time. It is provided to illustrate the features of the extensions, not to illustrate why you might need them. Much larger instances are are required for that.
1026
1027\section{Solving sub-problems in Parallel \label{parallelsec}}
1028
1029Pyomo makes use of the pyro facility to enable sub-problems to easily be assigned to be solved in parallel. This capability is suppored by pysp. We will refer to a single master computer
1030and multiple slave computers in the discussion here, but actually, the master computing processes can (and for synchronous parallel, should) be on a
1031processor that also runs a slave process.
1032
1033Here are the commands in order:
1034\begin{enumerate}
1035\item On the master: \verb|pyro-ns|
1036\item On the master: \verb|dispatch_srvr|
1037\item On each slave: \verb|pyro_mip_server|
1038\item On the master: \verb|runph ... --solver-manager=pyro ...|
1039\end{enumerate}
1040Note that the command \verb|pryo-ns| and the argument \verb|solver-manger| have a dash in the middle, while
1041the commands \verb|dispatch_srvr| and \verb|pyro_mip_server| have underscores. The first three commands launch processes
1042that have no internal mechanism for termination; i.e., they will be terminated only if they crash or if they are
1043killed by an external process. It is common to launch these processes with output redirection, such as \verb|pyro-ns >& pyrons.log|. The
1044\verb|runph| command is a normal runph command with the usual arguments with the additional specification that subproblem solves should
1045be directed to the pyro solver manager.
Note: See TracBrowser for help on using the repository browser.