[DEV] Define final cost in SPModel

Final cost is now defined in SPMOdel as
a PolyhedralFunction

Author: frapac

src/SDDPoptimize.jl

@@ -40,26 +40,18 @@ Returns :
 function solve_SDDP(model::SPModel,
                     param::SDDPparameters,
                     display=0::Int64,
-                    V=nothing,
-                    V_final=nothing)
+                    V=nothing)
 
-    # First step: process terminal costs.
-    # If not specified, default value is null functions
-    if isa(V_final, Void)
-        Vf = PolyhedralFunction(zeros(1), zeros(1, model.dimStates), 1)
-    else
-        Vf = V_final
-    end
-
-    # Second step: process value functions if hotstart is called
+    # First step: process value functions if hotstart is called
     if isa(V, Vector{PolyhedralFunction})
         # If V is already specified, then call hotstart:
         problems = hotstart_SDDP(model, param, V)
     else
         # Otherwise, initialize value functions:
-        V, problems = initialize_value_functions(model, param, Vf)
+        V, problems = initialize_value_functions(model, param)
     end
 
+    # Run SDDP upon example:
     run_SDDP!(model, param, V, problems, display)
     return V, problems
 end
@@ -207,8 +199,6 @@ Parameter:
 
 """
 function build_terminal_cost!(model::SPModel, problem::JuMP.Model, Vt::PolyhedralFunction)
-    alpha = getVar(problem, :alpha)
-
     # if shape is PolyhedralFunction, build terminal cost with it:
     alpha = getVar(problem, :alpha)
     x = getVar(problem, :x)
@@ -307,12 +297,9 @@ Return:
     the initialization of linear problems used to approximate
     each value function
 
-- V_final (PolydrehalFunction)
-    Terminal cost to penalize final state.
-
 """
 function initialize_value_functions( model::SPModel,
-                                     param::SDDPparameters, V_final::PolyhedralFunction)
+                                     param::SDDPparameters)
 
     solverProblems = build_models(model, param)
     solverProblems_null = build_models(model, param)
@@ -323,7 +310,10 @@ function initialize_value_functions( model::SPModel,
     # Build scenarios according to distribution laws:
     aleas = simulate_scenarios(model.noises, param.forwardPassNumber)
 
-    V[end] = V_final
+    V[end] = model.finalCost
+
+    # Add final costs to solverProblems:
+    build_terminal_cost!(model, solverProblems[end], V[end])
 
     stockTrajectories = forward_simulations(model,
                         param,
@@ -332,8 +322,6 @@ function initialize_value_functions( model::SPModel,
                         aleas,
                         false, true, false)[2]
 
-    build_terminal_cost!(model, solverProblems[end], V[end])
-
     backward_pass!(model,
                   param,
                   V,

src/forwardBackwardIterations.jl

@@ -111,6 +111,9 @@ function forward_simulations(model::SPModel,
 
             if returnCosts
                 costs[k] += nextstep.cost - nextstep.cost_to_go
+                if t==T-1
+                    costs[k] += nextstep.cost_to_go
+                end
             end
         end
     end

src/objects.jl

@@ -11,6 +11,15 @@ include("noises.jl")
 abstract SPModel
 
 
+type PolyhedralFunction
+    #function defined by max_k betas[k] + lambdas[k,:]*x
+    betas::Vector{Float64}
+    lambdas::Array{Float64,2} #lambdas[k,:] is the subgradient
+    # number of cuts:
+    numCuts::Int64
+end
+
+
 type LinearDynamicLinearCostSPmodel <: SPModel
     # problem dimension
     stageNumber::Int64
@@ -28,18 +37,28 @@ type LinearDynamicLinearCostSPmodel <: SPModel
     dynamics::Function
     noises::Vector{NoiseLaw}
 
-    function LinearDynamicLinearCostSPmodel(nstage, ubounds, x0, cost, dynamic, aleas)
+    finalCost::PolyhedralFunction
+
+    function LinearDynamicLinearCostSPmodel(nstage, ubounds, x0, cost, dynamic, aleas, Vfinal=nothing)
 
         dimStates = length(x0)
         dimControls = length(ubounds)
         dimNoises = length(aleas[1].support[:, 1])
 
+        # First step: process terminal costs.
+        # If not specified, default value is null function
+        if isa(Vfinal, Void)
+            Vf = PolyhedralFunction(zeros(1), zeros(1, dimStates), 1)
+        else
+            Vf = Vfinal
+        end
+
         xbounds = []
         for i = 1:dimStates
             push!(xbounds, (-Inf, Inf))
         end
 
-        return new(nstage, dimControls, dimStates, dimNoises, xbounds, ubounds, x0, cost, dynamic, aleas)
+        return new(nstage, dimControls, dimStates, dimNoises, xbounds, ubounds, x0, cost, dynamic, aleas, Vf)
     end
 end
 
@@ -60,17 +79,24 @@ type PiecewiseLinearCostSPmodel <: SPModel
     costFunctions::Vector{Function}
     dynamics::Function
     noises::Vector{NoiseLaw}
+    finalCost::PolyhedralFunction
 
-    function PiecewiseLinearCostSPmodel(nstage, ubounds, x0, costs, dynamic, aleas)
+    function PiecewiseLinearCostSPmodel(nstage, ubounds, x0, costs, dynamic, aleas, Vfinal=nothing)
         dimStates = length(x0)
         dimControls = length(ubounds)
         dimNoises = length(aleas[1].support[:, 1])
 
+        if isa(Vfinal, Void)
+            Vf = PolyhedralFunction(zeros(1), zeros(1, dimStates), 1)
+        else
+            Vf = Vfinal
+        end
+
         xbounds = []
         for i = 1:dimStates
             push!(xbounds, (-Inf, Inf))
         end
-        return new(nstage, dimControls, dimStates, dimNoises, xbounds, ubounds, x0, costs, dynamic, aleas)
+        return new(nstage, dimControls, dimStates, dimNoises, xbounds, ubounds, x0, costs, dynamic, aleas, Vf)
     end
 end
 
@@ -196,13 +222,6 @@ function set_max_iterations(param::SDDPparameters, n_iter::Int)
 end
 
 
-type PolyhedralFunction
-    #function defined by max_k betas[k] + lambdas[k,:]*x
-    betas::Vector{Float64}
-    lambdas::Array{Float64,2} #lambdas[k,:] is the subgradient
-    # number of cuts:
-    numCuts::Int64
-end
 
 
 type NextStep

test/runtests.jl

@@ -123,6 +123,8 @@ facts("SDDP algorithm: 1D case") do
         V, pbs = solve_SDDP(model, params, 0)
         @fact typeof(V) --> Vector{StochDynamicProgramming.PolyhedralFunction}
         @fact typeof(pbs) --> Vector{JuMP.Model}
+        @fact length(pbs) --> n_stages - 1
+        @fact length(V) --> n_stages
 
         # Test if the first subgradient has the same dimension as state:
         @fact length(V[1].lambdas[1, :]) --> model.dimStates