Some user-selectable engineering methods require a mixed integer programming solution. The Unit Power Table method, for example, requires integer programming to model the start-up and shut-down of individual turbines. Mixed integer programming adds substantial computational time to an Optimization model.

Inquiries can be sent to riverware-support@colorado.edu.

Mixed Integer Programs are optimization problems that contain both continuous and integer variables. In the river and reservoir modeling context, many of the variables are continuous. For example, outflow, pool elevation, and storage are all continuous variables. In contrast, hydropower turbines are operated such that a unit is either off or fully on. Thus, the Unit Power Table method has introduced some integer and binary (zero or one) variables. More specifically, the number of units generating at a given plant during a given time step is an integer variable and the decision to operate a single unit is a binary variable; the unit is either off or on. The modeling of integer and binary variables is sufficiently similar that the term integer is used to refer to both integer and binary variables.

Integer variables are also used to model discontinuities. For example, avoidance and cavitation zones are regions where power cannot be generated without damaging the equipment. These zones create a discontinuity in the power curve, and integer variables are used to model this discontinuity.

RiverWare solves this using Mixed Integer Linear Programs, but we will generally omit the word linear. The linear restriction reflects the fact that all non-linear functions have been replaced by linear or piecewise linear functions. The use of linear approximation allows the optimization to solve faster. We will also omit the word mixed and refer to Mixed-integer Linear Programming as simply Integer Programming.

Integer Programs are more difficult to solve than linear programs. For this reason, RiverWare has limited use of integer variables and the heuristic solution methods used are discussed in more detail in the rest of this section.

Two conditions are necessary for Integer Programming in RiverWare; integer programming is only allowed in the lowest priority policy and the policy must directly or indirectly use integer variables.

Integer Programming is limited to the lowest priority for performance reasons. To include integer variables at additional priorities would lead to what it called Integer Goal Programming (IGP). While there is no theoretical barrier, IGP is too computational challenging at this time. Consequently, RiverWare only allows integer variables in the final priority. If integer variables are referenced within a previous objective, RiverWare will relax the integrality constraints and solve the relaxed linear programming problem.

If a RiverWare policy doesn’t reference integer variables then integer programming is not used. Currently the only source of integer variables is the Unit Power Table and related methods (see “Unit Power Table”). Therefore, if only plant power methods are selected, integer variables will not be introduced and RiverWare will be solved as a Linear Goal Program.

The integer variables can be introduced either by directly referencing them in policy or by indirectly referencing them in policy. Integer variables are typically introduced indirectly by a user policy that optimizes the economic value of power using slots on the model’s Thermal object. During optimization, these slot variables will be replaced with equivalent expressions involving Power and related quantities on the individual power reservoirs. In this way, if the Unit Power Table method is selected, constraints which model unit level power concerns will be introduced into the optimization problem.

In addition, the user might write policy constraints which refer to binary quantities. For example, the user might want to forbid scheduling holes; that is, turning a unit off for a single time period. One normally thinks of writing this constraint as

However, to avoid referencing future time periods they will need to write this as:

Understanding the method used to solve integer programs is useful for understanding the solution information and adjusting parameters that control the solution process.

Integer programs are solved by first starting with a linear programming relaxation, for which integrality restrictions (i.e. that variables are either 0 or 1) are temporarily removed. From this relaxed solution two subproblems are created by branching on an integer variable. For example, the solver might branch on Unit 1 being used in the first period, with one branch requiring it to be used and the other branch requiring it not be used. This branching process can be continued with additional variables. The resulting process creates a tree structure that contains an optimal solution at one of the leaf nodes. Explicitly exploring the whole tree is impossible for all but the smallest problems. Fortunately, large portions of the tree can frequently be pruned off without exploration. The pruning is possible because the linear programming relaxation of a node in the tree may either be an integer solution, infeasible, or sufficiently suboptimal.

The integer programs formulated in RiverWare are sufficiently hard that an optional, but highly recommended heuristic has been added to speed up the solution process. The heuristic is an iterative timestep approach that focuses on solving early timesteps with integer variables while temporarily relaxing the integrality conditions for later timesteps. By relaxing the integrality constraints for most of the timesteps, the problem can be solved much more quickly. After obtaining a solution, the integer variables for the earliest timesteps are frozen at their integer values, and some previously relaxed variables are treated as integer variables. The problem is resolved, and the process continues until all timesteps have been modeled as integer variables and frozen. For each iteration the integer variables are constrained both by the frozen integer variables for previous time steps and the linear programming relaxation for future timesteps.

There are two time windows that are advanced at each iteration in this process, the range of (non-frozen) integer variables in the problem, and the possibly smaller range of integer variables that should be frozen after the problem is solved. The size of each of these windows is an important parameter of this process. By default, one time step is used for each of these time windows, but the user can increase these values by setting appropriate parameters, Mip Integer Window Size and Mip Freeze Window Size. If both variables are set to the number of timesteps in the run, then the MIP will be solved as a single problem containing all of the integer variables. At this time, we do not recommend increasing the window sizes beyond 1 because of the dramatic increase in computation time.

This process may find the optimal solution, but it is not guaranteed. During testing, the timestep approach has found high quality, feasible solutions in a fraction of the time to find an optimal solution and prove it is optimal.

See “Mathematical Formulation of MIP for Unit Power” for details on the mathematical formulation of the mixed integer programming to model unit power.

CPLEX provides a variety of parameters to tune the solution of optimization problems. RiverWare makes these parameters accessible to the user. In addition, RiverWare has added some additional parameters to the list that are related to RiverWare’s use of goal programming and time windows. All of the parameters have default values determined by CPLEX or by CADSWES. In some cases CADSWES has replaced the default values of CPLEX parameters with values that have performed better for RiverWare models. These changes in the CPLEX parameters are in the cplex.par file packaged with RiverWare. Similarly, any non-default values for CADSWES parameters are contained in the goal.par file.

The vast majority of the parameters need never be examined, but some of them are potentially useful for integer programming. While branch and bound will theoretically find an optimal solution given enough time, the exponential growth of the search tree means that the time to do so is prohibitive for many practical problems. Instead, the tree can be used in a heuristic manner to find near optimal solutions. By limiting the number of solutions, the time allowed, and/or the optimality tolerance, the search can be terminated early with a good solution.

The parameters (both for CPLEX and RiverWare specific parameters) are set through the RiverWare interface. From the Run Control dialog, select Optimization as the controller and then select the menu option View > Optimization Run Parameters. From the Optimization Run Parameters dialog, select Set CPLEX and Goal Parameters. The parameters are organized by categories. Select the desired parameter and change the current value to the desired value. Select either Apply (to save and continue with other parameter changes) or Ok (to save and exit the parameter dialog). If you check the Save non-default settings to parameter file box these settings will be saved to a file called cplex.init.par, and this file will be used to initialize parameter settings for future RiverWare sessions.

This parameter usually will not be effective when the timestep approach is used (see “Timestep Approach”). However, when optimizing the entire problem, this parameter can be very influential because of the structure of the branch and bound solution method. For maximization problems, the lower cutoff parameter will help prune the search tree until an integer feasible solution is found (an upper cutoff parameter can be used for minimization problems). Ideally, this parameter will be set high enough to eliminate large portions of the tree without eliminating all of the feasible integer solutions. Since the optimal integer solution isn’t known in advance, use of this parameter requires a little experimentation.

A starting point for the experimentation is the linear programming relaxation. By definition, integer solutions cannot improve on this value. For a maximization problem this means the optimal integer solution will be less than or equal to the linear programming relaxation and the lower cutoff parameter should certainly be set less than the objective function value for the linear programming relaxation. Our limited experimentation with test models suggests cutoff values ranging from about 90% of the linear programming relaxation to just below the linear programming relaxation.

Some experimentation is unavoidable with this parameter. If the lower cutoff value is set too high the search will end without a feasible solution. If the value is set too low (or there are too many time steps) the run will continue until some parameter stops the run or the user terminates the run. In both cases there won’t be an integer solution. If the problem is small enough to be solved, a reasonable cutoff value can typically be determined with a couple of tries.

The objective function of the linear programming relaxation can be obtained approximately by terminating the run with a little delay after the last objective has started. The delay is necessary to give CPLEX time to solve the relaxation. After terminating, the objective function of the best node in the tree will be printed in the RiverWare diagnostics window. This may be slightly less than the linear programming relaxation. In general this value will decrease as branch and bound proceeds. The number of nodes processed in the tree are also printed in the diagnostic. Unfortunately, for technical reasons the value can’t be printed while CPLEX is still running. The objective function is printed in optimization units, which are typically scaled from user units for numerical stability. For this reason, direct interpretation of the value is not advised. However, these are the same units that should be used for the lower cutoff value.

The lower cutoff parameter is in the MIP Tolerance Parameters category.

Sometimes, branch and bound will find a good or even optimal solution early in an optimization and spend the bulk of its time proving optimality or making relatively minor improvements in the solution. Limiting the number of solutions can improve the solution time. However, testing has shown that setting this parameter overly low while using the timestep approach can lead to integer infeasibility in the middle of the timestep iterations.

The solutions parameter is in the MIP Limit Parameters category.

Optimization problems can take a long time to solve. A time limit parameter allows a solver to exit gracefully when the solver takes too long. For integer programs, exiting when the time limit expires will not generate an error if a feasible solution has been found. Thus, a time limit parameter can be used to truncate exploration of the branch and bound tree.

The time limit parameter in the General Parameters category. The default value is set to 3000 seconds (50 minutes).

Optimality tolerance parameters define how close the solution must be to optimal to be considered close enough, either as a fraction of optimality (0.0 - 1.0) or as an absolute number (in optimization units). Even on relatively easy problems, proving 100% optimality (an optimality tolerance of zero) can require very large computation times. A small optimality tolerance (e.g., e-4 to e-6 absolute) generally yields a great improvement in runtime, and for certain problems larger tolerances (e.g., 1-10% fractional) can also be beneficial. Experimenting with these parameters may prove worthwhile.

The fractional and absolute optimality tolerance parameters are in the MIP Tolerance Parameters category and are called respectively, mipgap and absmipgap.

At the end of optimization, the optimized values are known but haven’t been returned to the workspace yet. The challenge is to know what variables will be returned to drive the post-optimization RBS run. When using the Unit Power Table method, the Rulebased Simulation can be driven in several ways:

Optimization could return the individual Unit Turbine Releases and some other piece of information to cause the object to dispatch and solve. For example, the additional information could be the Outflow, plant Turbine Release or the “U” flag (implementation soon) on Turbine Release. In fact, it is recommended, to return the individual Unit Turbine Release and total Outflow. This approach calculates power very similarly to the optimization calculation and is perhaps the most natural way to do a post-optimization simulation. Also, given unit turbine releases will always produce an answer.

Optimization could return the individual unit energy, the total Spill, and some other piece of information to cause the object to dispatch and solve. For example, the additional information could be the Plant Energy or the “U” flag (implementation soon) on Energy. This approach would be similar to the first option, but would return unit energies instead of unit flows. In both of these options, given the individual unit values (and even spill) would not be enough to cause the reservoir to dispatch. A plant level value must be specified. Given unit energies, there is the potential to have an infeasible solution in simulation. An operating head and flow could be calculated that because of linearization errors, cannot meet the given energies. This approach is not supported or recommended.

For the time being, the first option is the only supported. The rule should return the plant level Outflow and the Unit Turbine Release, marking both with rule flags. The plant level value will cause dispatching in simulation. With these two approaches, the user has options on what should be returned from optimization and how the object should solve. Experimentation will allow the user to determine the variables that provides the best solution for the application.

This section presents an improved mathematical formulation of the relationship between power and other quantities. This formulation models power generation at the unit level to capture the fact that power generation occurs at one power reservoir with a heterogeneous collection of turbines. Integer variables allow the user to model the fact that these units have discrete states (off, generating power, spinning, etc.) and transitions between states are potentially costly.

Table A.1 summarizes information concerning the slots used by the Unit Power Method, methods in the related dependent categories, and methods on the thermal object. The “Symbol” column indicates the symbolic notation that will be used in subsequent sections to refer to the slot.

The slot dimensionality column is included to give a rough idea of the resources required by these methods. This column uses the following abbreviations:

At the heart of the optimization controller’s model of power is the triangulation algorithm for partitioning the power curve domain (specifically flow and head), which are described here. The method successfully excludes zones that are infeasible because of vibration, cavitation, etc.

The input for this method is a three-dimensional power table and a table indicating zones that must be avoided. It is assumed that the power table contains the end points of head that define the zone and for each head the flow values that define the zone. The points themselves are assumed to be feasible operating points unless they represent minimum or maximum values of head or flow. There may be multiple separate avoidance zones.

Figure A.1 is an example of a triangulated power domain with an avoidance zone. In this figure, each vertex corresponds the data in a Unit Power Table row which apply to a single unit, head (operating head or net head), and Turbine Release. Power, an additional dimension represented in that table, is not shown.

Let us identify the avoidance zone values from the power table in the following way for the algorithm:

The basic idea is to iterate through the Head values in the (reduced) Unit Power Table, partitioning the vertical space between that head and the next into triangles. In general when deciding which of two potential triangles to create, RiverWare creates the one that has the lesser slope for the upper side.

If head is specified, then the triangles become line segments, and line segments become points.

In the context of optimization, a lambda set is (roughly) a discretization of a continuous space, so RiverWare uses this technique to model the discretization of the power curve domain into triangles. While solvers don’t require it, usually it is assumed that lambda variables are constrained as follows:

Lambda sets may additionally be constrained to be a specially ordered set of type 1 or 2, aka SOS1 or SOS2. The individual variables in an SOS are ordered and a reference row of weights is associated with each variable. An SOS1 may have at most one value that is nonzero while an SOS 2 may have at most two nonzero values and they must be adjacent. Solvers internally treat these sets as being similar to integer variables with branching to enforce the set conditions. The branching is based on weights in the reference row. The general idea is to calculate a weighted average of the variables using the reference row. The two branches are roughly that either variables on one side of the weighted average are zero or the variables on the other side are zero.

In addition, frequently the individual variables belonging to an SOS1are also binary variables. The variables in a SOS2 are almost always continuous.

Power is a nonlinear function of both flow and head. Hence, when it is discretized, the domain of the function is two-dimensional. The SOS technique can be applied here as well. In this case, a combination of 3 points is used to define flow and head. Alexander Martin has described how to model this case as an SOS type k, where k is 3 in this case (A. Martin, “Approximation of Non-linear functions in Mixed Integer Programming” Workshop on Integer Programming and Continuous Optimization, Chemnitz University of Technology, November 7-9, 2004, www-user.tu-chemnitz.de/~helmberg/workshop04/martin.ppt). Martin then modifies a solver to directly branch on the SOS k, creating separate branches for each region of the domain.

A variant of his method that is compatible with CPLEX and other solvers is used here. The flow and head domains are partitioned into triangles and define a variable for each triangle, Lambda_1. The triangles in an SOS 3 are analogous to the line segments in an SOS 2. The triangulation is described in section 2.7.

The following lambda variables are part of the formulation:

This variable binary, technically an SOS k, Lambda_2 is used to implement this functionality. Therefore, there is no need for a reference row and its weights. This variable may represent 2 points instead of 3 if a feasible region is one-dimensional. For example most data sets will contain points with zero flow for several heads, but disallow small flows. Thus, the segment of zero flow between adjacent heads will be a Lambda_1 variable even though it is not a triangle. In theory, there could be a single isolated operating point that would need its own Lambda_1 variable, but this is unlikely in practice.

This section presents the constraints which will be introduced into the optimization problem as power-related quantities are encountered in the optimization policy (when the Unit Power method is selected).

All variables are constrained to be positive unless otherwise noted.

These variables are required to be zero or one, and will be enforced on user inputs in the corresponding slots at the beginning of the run.

The following symbols will be used in the formulation below (see also slot summary table):

The following indexes are used:

Sigma indicates capital sigma used for a summation. The arguments following sigma in parenthesis are respectively the indexes to sum over and the expression summed.

Bold typeface for a variable indicates the constraint is a defining constraint for that variable.

The constraints are presented with an unspecified timestep (t) and unit (u). Instantiations will be created during the run only for those units and timesteps which are referenced by the user’s policy.

Any plant variable, PVx, where x represents the appropriate subscripts is defined by the sum of the unit variables, Vux, indexed by unit u and x. Plant variables are only created if they are used directly or indirectly by policy. Thus, the defining constraint is as follows:

The modeling of power will be more accurate when lambda variables and SOS sets are used to replace the plant level linearization of pool elevation and tailwater. Both of these affect operating head. This has been temporarily postponed until a potential problem is addressed: the associated lambda variables might be frozen by an objective solved using linear programming before it reaches the power objective. If this were to occur, the intended model for integer programming would be disrupted. This is more likely to be a problem for pool elevation than tailwater because pool elevation is important for non-power policy, whereas tailwater is usually only important for power. Moreover, the error in existing methods is larger for tailwater than pool elevation. For these reasons, it is more likely beneficial to change the modeling of tailwater.

The pool elevation constraints are:

If the tailwater method is optTWStageFlowLookupTable, the tailwater constraints define tailwater, outflow, and tailwater base value in terms of Lambda_8 points directly or indirectly in terms of Lambda_9 and Lambda_10:

If the tailwater method is optTWBaseValuePlusLookupTable, the new tailwater constraints define temp tailwater lookup and outflow in terms of Lambda_11:

For this method, the existing tailwater constraint is retained:

The thermal object has a number of multi-slots that are potential variables in an optimization problem. These multi-slots typically represent the system wide total of individual reservoir values that are linked to these slots. Equations for these sums are added automatically if a policy either directly or indirectly references one of the thermal object multi-slots.

Many of the existing variables and definitions are already on the thermal object. In addition, to the existing equations, the following new equations will be added.

System Regulation Up is greater than the sum of Block Regulation Up variables:

The Regulation Up Value equals the sum of the value of the Block Regulation Up variables:

System Regulation Down is greater than the sum of Block Regulation Down variables:

The Regulation Down Value equals the sum of the value of the Block Regulation Down variables: