####################################
# New Pyomo DAE Components and Use #
####################################

Components
============
* DifferentialSet
* Differential

Discretization Transformations
===============================
* Implicit_Euler
* Collocation_Discretization

Component Descriptions and Use
================================

-------------------------------- 
DifferentialSet
-------------------------------- 

This component represents the term in the denominator of a derivative,
i.e. what the derivative is being taken with respect to or the
independent variable(In many cases this is time). Minimally, this set
must contain two indices indicating the initial and final values for
the independent variable where the initial value must be less than the
final value. Additional indicies may be specified if the user wants to
specify certain points in the range where finite element boundaries
must be located during the discretization. All indicies must be
numeric values. Advanced users may initialize the DifferentialSet with
all the desired discretization points before discretizing the model in
order to save processing time during the discretization
transformation.

--------------------------------

Example Use:
model.t = DifferentialSet(bounds=(0,5))

model.t = DifferentialSet(initialize=[0,1,2,5])

model.t = DifferentialSet()
# .dat file
set t := 0 0.5 2.25 3.75 5;

** Notice that if a .dat file is used to initialize a DifferentialSet
it is done using the 'set' command and not 'differentialset'**

---------------------------------
Differential
---------------------------------

This is a component used to specify differential equations in the
model. The user must write a function which returns the right hand
side of an equation in the form 
dx/dt = f(x,y,u,t,p) 
where
x: differential variables
y: algebraic variables
u: control variables
t: independent variable
p: parameters which do not depend on the independent variable

A differential has no positional arguments. The required keyword
arguments are 'dvar' and 'rule'. 'dvar' represents the differential
variable or the variable whos derivative is equal to the expression
specified in the rule. The differential variable must be *explicitly*
indexed by at least one DifferentialSet. If the differential variable
is indexed by multiple DifferentialSets then the 'dset' keyword
argument must be specified in the Differential to indicate which
DifferentialSet is to be used as the independent variable in the
differential equation. If the differential is only indexed by a single
DifferentialSet then that DifferentialSet is assumed to be the
independent variable and the 'dset' keyword argument is optional.

The 'rule' keyword argument should be set to a function which returns
the expression for the right hand side of the differential
equation. The indexing on the function should exactly match the
indexing on the differential variable in order and
number. Additionally, the indexing on the Differential object will
match that of the differential variable.

An optional 'bounds' keyword argument may be specified if one would
like to set bounds on the derivative of the differential
variable. 'bounds' can be set to a tuple if they are constant for all
derivatives or they can be set to a function which computes the bounds
for a particular derivative.

An optional 'initialize' keyword argument may be specified if one
would like to set the initial value of the derivative. Note, this is
not the same as setting the initial condition. This functionality is
provided since the starting point can play a large part in being able
to solve a nonlinear programming problem. 'initialize' can be set to a
constant, a function, or a dictionary.

--------------------------------

Example Use:
model.t = DiffSet(bounds=(0,5))
model.tray = RangeSet(8)
model.x1 = Var(model.t)
model.x2 = Var(model.tray, model.t)

def _x1dot(model, ti):
    return model.x1[ti]**2 + ti*6
model.x1dot = Differential(dset=model.t,dvar=model.x1, rule=_x1dot)

def _x2dot(model, i, ti):
    return model.x2[i,ti]*3 + ti**2
model.x2dot = Differential(dvar=model.x2, rule=_x2dot, bounds=(None,5),
	    initialize=3)

If you would like to use a derivative in a separate constraint (such
as setting an initial or final condition) you can access it like a regular
variable. See the example below.

Example Accessing a Derivative (Using the Differentials previously defined):

def _x1dot_final(model):
    return model.x1dot[5] == 0
model.x1dot_final = Constraint(rule=_x1dot_final)

def _x2dot_init(model,i):
    return model.x2dot[i,0] == 0
model.x2dot_init = Constraint(model.tray,rule=_x2dot_init)

======================================================
Discretization Transformations: Descriptions and Use
======================================================

Before a Pyomo model with Differential components can be sent to a
solver it must first be sent through a discretization
transformation. These transformations convert any differential
equations in the model into nonlinear equations by using various
discretization techniques. This is known as the simultaneous approach
for numerically solving systems of differential algebraic
equations. Two discretization schemes have been implemented in Pyomo,
Implicit Euler and Orthogonal Collocation. The user must write a
Python script in order to use these transformations. Examples of the
required scripting steps are shown below for each of the
discretization schemes.

A transformation framework along with certain utility functions have
been created so that advanced users may easily implement custom
discretization schemes other than the two already implemented in
Pyomo. The tranformation framework consists of the following steps:

1) Set Discretization Options
2) Discretize the Differential Set(s)
3) Update Model Components
4) Add Discretization Equations
5) Return Discretized Model

The transformations return valid Pyomo models which can be further
manipulated before being sent to a solver. Examples of this are shown
in the 'Example Use' sections of both transformations.

---------------------------------- 
Implicit Euler
---------------------------------- 

This transformation uses the Implicit Euler (also called Backward
Euler) method to discretize the differential equations in the
model. The discretization equations associated with this method are
shown below:

dx/dt = f(t,x)
discretize t and x such that
x(t0) = x_{0}
x(t0+kh) = x_{k}
x_{k+1} = x_{k}+h*f(t_{k+1},x_{k+1})
t_{k+1} = t_{k}+h

where h is the step size between discretization points or the size of
each finite element. These are equations/constraints that are
generated automatically when the Implicit_Euler transformation is
applied to a Pyomo model. 

There are two discretization options available to the Implicit_Euler
transformation which can be specified as keyword arguments to the
'apply' function of the transformation object. The first is
'nfe'. This keyword argument is used to specify the number of finite
elements to be included in the discretization, i.e. it gives the user
some indirect control over the value of h. The default value for 'nfe'
is 10. If the existing number of finite elements in a differential set
is less than the desired number, new discretization points will be
added to the set. *NOTE: Discretization points will not be removed
from the set if the existing number of finite elements exceeds the
number specified in the keyword argument*

The other keyword argument for the ImplicitEuler transformation is
'diffset'. This is an optional keyword and allows the user to specify
which DifferentialSet is to be discretized. This is for the case when
the user has multiple DifferentialSets in a model which should be
discretized differently. If the 'diffset' keyword argument is not used
then the same discretization scheme will be applied to every
DifferentialSet in the model.

--------------------------------

Example Use: 
The following code is a Python script applying the Implicit_Euler
discretization. The code also shows how to add a constraint to a
discretized model.

from coopr.pyomo import *
from coopr.dae import *
from coopr.opt import SolverFactory
from coopr.dae.impliciteuler import Implicit_Euler

# Import Pyomo model
from pyomoExample import model

# Create model
instance = model.create()

# Discretize model using Implicit Euler method
discretize = Implicit_Euler(instance)
disc_instance = discretize.apply(nfe=20,diffset=instance.time)

# Add another constraint to discretized model
def _sum_limit(model):
    return sum(model.x1[i] for i in model.time) <= 50
disc_instance.con_sum_limit = Constraint(rule=_sum_limit)

# Solve discretized model
opt = SolverFactory('ipopt')
results = opt.solve(disc_instance)

# Load Results
disc_instance.load(results)
 
----------------------------------
Collocation
----------------------------------

This transformation uses orthogonal collocation to discretize the
differential equations in the model. Currently, Gauss-Radau
collocation is the only type that has been implemented. For more
information on orthogonal collocation and the discretization equations
associated with this method please see chapter 10 of the book
"Nonlinear Programming: Concepts, Algorithms, and Applications to
Chemical Processes" by L.T. Biegler.

There are three discretization options available to the
Collocation_Discretization transformation which can be specified as
keyword arguments to the 'apply' function of the transformation
object. The first two are 'nfe' and 'diffset'. Please see the
descriptions of these options in the section on the Implicit_Euler
transformation as they work the same in this case. The last option
available to the collocation transformation is 'ncp' which can be used
to specify the number of collocation points within each finite
element. The default value is 3. If the user's version of Python has
access to the package Numpy then any number of collocation points may
be specified, otherwise the maximum number is 10. 

*NOTE: Any points that exist in a DifferentialSet before
 discretization will be used as finite element boundaries and not as
 collocation points. The locations of the collocation points cannot be
 specified by the user, they must be generated by the transformation*

In some applications it may be desirable to use a different number of
collocation points for certain variables in the model. For example, in
optimal control problems, the control variable is often restricted to
be constant over each finite element while the other state variables
are not. To formulate a problem such as this, the user can use the
'reduce_collocation_points' function. To use this function the
following keyword arguments must be specified: 

'var': The variable being restricted to fewer collocation points
'diffset': The differential set indexing the specified 'var' which has
been previously discretized.
'ncp': The new number of collocation points. Must be at least 1 and
less than the number of collocation points used to discretize the
'diffset'.

The 'reduce_collocation_points' function may only be applied to a
model after the 'apply' function has been called to discretize the
specified DifferentialSet. The function works by adding constraints
to the discretized model which force any extra, undesired collocation
points to be interpolated from the others.

--------------------------------

Example Use: 
The following code is a Python script applying the
Collocation_Discretization transformation. The code also shows how to
use the 'reduce_collocation_points' function and how to add an
objective function to a discretized model.

from coopr.pyomo import *
from coopr.dae import *
from coopr.opt import SolverFactory
from coopr.dae.colloc import Collocation_Discretization

# Import Pyomo model
from pyomoExample2 import model

# Create model
instance = model.create()

# Discretize model using Radau Collocation
discretize = Collocation_Discretization(instance)
disc_instance = discretize.apply(nfe=20, ncp=6)

# Control variable u made constant over each finite element
disc_instance = discretize.reduce_collocation_points(
    var=instance.u, ncp=1, diffset=instance.time)

# Add objective function after model has been discretized
def obj_rule(model):
    return sum((model.x[i]-model.x_ref)**2 for i in model.time)
disc_instance.obj = Objective(rule=obj_rule)

# Solve discretized model
opt = SolverFactory('ipopt')
results = opt.solve(disc_instance)

# Load Results
disc_instance.load(results)
