[WIP] Improving initialization

REQUIRE

@@ -1,5 +1,7 @@
 julia 0.2-
+Compat
 Docile
+JuMP
 Reexport
 Requires
 Sundials

examples/basics/vanderpol_with_events.jl

@@ -12,7 +12,7 @@ An extension of Sims.Examples.Basics.Vanderpol. Events are triggered
 every 2 sec that change the quantity `mu`.
 """ ->
 function VanderpolWithEvents()
-    y = Unknown(1.0, "y")   
+    y = Unknown(value = 1.0, label = "y", fixed = true)   
     x = Unknown("x")
     mu = Discrete(1.0)
     alpha = @liftd(0.8 * :mu)
@@ -21,8 +21,7 @@ function VanderpolWithEvents()
     alpha_u = Unknown(value(alpha), "alpha_u") 
     beta_u = Unknown(value(beta), "beta_u") 
     @equations begin
-        # The -1.0 in der(x, -1.0) is the initial value for the derivative 
-        der(x, -1.0) = mu * (1 - y^2) * x - y
+        der(x) = mu * (1 - y^2) * x - y
         der(y) = x
         mu_u    = mu
         alpha_u = alpha

examples/lib/Lib.jl

@@ -42,7 +42,7 @@ function runexamples()
 
     ## Power systems
     rlm   = sim(RLModel(), 0.2)
-    pim   = sim(PiModel(), 0.02)
+    ## pim   = sim(PiModel(), 0.02)
     ## mm    = sim(ModalModel(), 0.2)
 
     ## Rotational

examples/lib/electrical.jl

@@ -272,7 +272,7 @@ voltages).
 function ChuaCircuit()
     n1 = Voltage("n1")
     n2 = Voltage("n2")
-    n3 = Voltage(4.0, "n3")
+    n3 = Voltage(value = 4.0, label = "n3", fixed = true)
     g = 0.0
     function NonlinearResistor(n1::ElectricalNode, n2::ElectricalNode, Ga, Gb, Ve)
         i = Current(compatible_values(n1, n2))

src/Sims.jl

@@ -16,7 +16,7 @@ import Base.ifelse,
 ## Types
 export ModelType, UnknownCategory, Unknown, UnknownVariable, DefaultUnknown, DerUnknown, RefUnknown, 
        RefBranch, InitialEquation, Model, MExpr, Event, LeftVar, StructuralEvent,
-       EquationSet, SimFunctions, Sim, SimResult
+       EquationSet, SimFunctions, Sim, SimState, SimResult
 
 export UnknownReactive, Discrete, Parameter
 
@@ -27,7 +27,8 @@ export Equation, @equations, is_unknown, der, delay, mexpr, compatible_values, r
        basetypeof, from_real, to_real,
        gplot, wplot,
        check, sim_verbose, 
-       elaborate, create_sim, create_simstate, sim, sunsim, dasslsim, solve
+       elaborate, create_sim, create_simstate, sim, sunsim, dasslsim, solve,
+       initialize!
 
 ## Model methods
 export Branch, BoolEvent
@@ -50,6 +51,7 @@ export DocTag
 ##     :( @doc $ex -> type $(gensym()) <: DocTag end )
 ## end
 ## export @doctag, DocTag
+using Compat
 
 
 include("main.jl")

src/dassl.jl

@@ -53,20 +53,24 @@ function dasslrootfun(neq, t_in, y_in, yp_in, nrt, rval_out, rpar, ipar)
     return nothing
 end
 
+initdassl = @compat Dict(:none => 0, :Ya_Ydp => 1, :Y => 2)
 
 @doc* """
 The solver that uses DASKR, a variant of DASSL.
 
 See [sim](#sim) for the interface.
 """ ->
-function dasslsim(ss::SimState, tstop::Float64=1.0, Nsteps::Int=500, reltol::Float64=1e-4, abstol::Float64=1e-4)
+function dasslsim(ss::SimState, tstop::Float64=1.0, Nsteps::Int=500, reltol::Float64=1e-4, abstol::Float64=1e-4, init::Symbol=:Ya_Ydp)
     # tstop & Nsteps should be in options
     sim_info("starting sim()")
 
     sm = ss.sm
     for x in sm.discrete_inputs
-        push!(x, x.initialvalue)
+        push!(x.signal, x.initialvalue)
     end
+    ss.y[:] = ss.y0
+    @show ss.y
+    ss.yp[:] = ss.yp0
     yidx = sm.outputs .!= ""
     ## yidx = map((s) -> s != "", sm.outputs)
     Noutputs = sum(yidx)
@@ -75,14 +79,14 @@ function dasslsim(ss::SimState, tstop::Float64=1.0, Nsteps::Int=500, reltol::Flo
     tout = [tstep]
     idid = [int32(0)]
     info = fill(int32(0), 20)
-    info[11] = 1    # calc initial conditions (1 or 2) / don't calc (0)
+    info[11] = initdassl[init]    # calc initial conditions (1 or 2) / don't calc (0)
     info[18] = 2    # more initialization info
 
     function setup_sim(ss::SimState, tstart::Float64, tstop::Float64, Nsteps::Int; reltol::Float64=1e-5, abstol::Float64=1e-3)
         N = [int32(length(ss.y0))]
         t = [tstart]
-        y = copy(ss.y0)
-        yp = copy(ss.yp0)
+        y = ss.y
+        yp = ss.yp
         sm = ss.sm
         nrt = [int32(length(sm.F.event_pos))]
         rpar = [0.0]
@@ -167,13 +171,13 @@ function dasslsim(ss::SimState, tstop::Float64=1.0, Nsteps::Int=500, reltol::Flo
 
                     # Restart the simulation:
                     info[1] = 0
-                    info[11] = 1    # do/don't calc initial conditions
+                    info[11] = initdassl[init]    # do/don't calc initial conditions
                     simulate = setup_sim(ss, t[1], tstop, int(Nsteps * (tstop - t[1]) / tstop), reltol=reltol, abstol=abstol)
                     yidx = sm.outputs .!= ""
                 elseif any(jroot .!= 0)
                     sim_info("event found at t = $(t[1]), restarting")
                     info[1] = 0
-                    info[11] = 1    # do/don't calc initial conditions
+                    info[11] = initdassl[init]    # do/don't calc initial conditions
                 end
             end
         elseif idid[1] < 0 && idid[1] > -11
@@ -186,20 +190,20 @@ function dasslsim(ss::SimState, tstop::Float64=1.0, Nsteps::Int=500, reltol::Flo
     end
     SimResult(yout, [sm.outputs[yidx]])
 end
-dasslsim(ss::SimState; tstop = 1.0, Nsteps = 500, reltol = 1e-4, abstol = 1e-4) =
-    dasslsim(ss, tstop, Nsteps, reltol, abstol)
+dasslsim(ss::SimState; tstop = 1.0, Nsteps = 500, reltol = 1e-4, abstol = 1e-4, init = :Ya_Ydp) =
+    dasslsim(ss, tstop, Nsteps, reltol, abstol, init)
 
-dasslsim(m::Model, tstop = 1.0, Nsteps = 500, reltol = 1e-4, abstol = 1e-4) =
-    dasslsim(create_simstate(m), tstop, Nsteps, reltol, abstol)
-dasslsim(m::Model; tstop = 1.0, Nsteps = 500, reltol = 1e-4, abstol = 1e-4) =
-    dasslsim(create_simstate(m), tstop, Nsteps, reltol, abstol)
+dasslsim(m::Model, tstop = 1.0, Nsteps = 500, reltol = 1e-4, abstol = 1e-4, init = :Ya_Ydp) =
+    dasslsim(create_simstate(m), tstop, Nsteps, reltol, abstol, init)
+dasslsim(m::Model; tstop = 1.0, Nsteps = 500, reltol = 1e-4, abstol = 1e-4, init = :Ya_Ydp) =
+    dasslsim(create_simstate(m), tstop, Nsteps, reltol, abstol, init)
 
-dasslsim(sm::Sim, tstop = 1.0, Nsteps = 500, reltol = 1e-4, abstol = 1e-4) =
-    dasslsim(create_simstate(sm), tstop, Nsteps, reltol, abstol)
-dasslsim(sm::Sim; tstop = 1.0, Nsteps = 500, reltol = 1e-4, abstol = 1e-4) =
-    dasslsim(create_simstate(sm), tstop, Nsteps, reltol, abstol)
-
-dasslsim(e::EquationSet, tstop = 1.0, Nsteps = 500, reltol = 1e-4, abstol = 1e-4) =
-    dasslsim(create_simstate(e), tstop, Nsteps, reltol, abstol)
-dasslsim(e::EquationSet; tstop = 1.0, Nsteps = 500, reltol = 1e-4, abstol = 1e-4) =
-    dasslsim(create_simstate(e), tstop, Nsteps, reltol, abstol)
+dasslsim(sm::Sim, tstop = 1.0, Nsteps = 500, reltol = 1e-4, abstol = 1e-4, init = :Ya_Ydp) =
+    dasslsim(create_simstate(sm), tstop, Nsteps, reltol, abstol, init)
+dasslsim(sm::Sim; tstop = 1.0, Nsteps = 500, reltol = 1e-4, abstol = 1e-4, init = :Ya_Ydp) =
+    dasslsim(create_simstate(sm), tstop, Nsteps, reltol, abstol, init)
+
+dasslsim(e::EquationSet, tstop = 1.0, Nsteps = 500, reltol = 1e-4, abstol = 1e-4, init = :Ya_Ydp) =
+    dasslsim(create_simstate(e), tstop, Nsteps, reltol, abstol, init)
+dasslsim(e::EquationSet; tstop = 1.0, Nsteps = 500, reltol = 1e-4, abstol = 1e-4, init = :Ya_Ydp) =
+    dasslsim(create_simstate(e), tstop, Nsteps, reltol, abstol, init)

src/main.jl

@@ -398,31 +398,31 @@ end
 
 # special case to avoid a warning:
 import Base.(^)
-(^)(x::ModelType, y::Integer) = mexpr(:call, (^), _expr(x), y)
+(^)(x::ModelType, y::Integer) = mexpr(:call, :(^), _expr(x), y)
 
 for f in binary_functions
     ## @eval import Base.(f)
     eval(Expr(:toplevel, Expr(:import, :Base, f)))
-    @eval ($f)(x::ModelType, y::ModelType) = mexpr(:call, ($f), _expr(x), _expr(y))
-    @eval ($f)(x::ModelType, y::Number) = mexpr(:call, ($f), _expr(x), y)
-    @eval ($f)(x::ModelType, y::AbstractArray) = mexpr(:call, ($f), _expr(x), y)
-    @eval ($f)(x::Number, y::ModelType) = mexpr(:call, ($f), x, _expr(y))
-    @eval ($f)(x::AbstractArray, y::ModelType) = mexpr(:call, ($f), x, _expr(y))
+    @eval ($f)(x::ModelType, y::ModelType) = mexpr(:call, $(Expr(:quote, f)), _expr(x), _expr(y))
+    @eval ($f)(x::ModelType, y::Number) = mexpr(:call, $(Expr(:quote, f)), _expr(x), y)
+    @eval ($f)(x::ModelType, y::AbstractArray) = mexpr(:call, $(Expr(:quote, f)), _expr(x), y)
+    @eval ($f)(x::Number, y::ModelType) = mexpr(:call, $(Expr(:quote, f)), x, _expr(y))
+    @eval ($f)(x::AbstractArray, y::ModelType) = mexpr(:call, $(Expr(:quote, f)), x, _expr(y))
 end 
 
 for f in unary_functions
     ## @eval import Base.(f)
     eval(Expr(:toplevel, Expr(:import, :Base, f)))
 
     # Define separate method to get rid of 'ambiguous definition' warnings in base/floatfuncs.jl
-    @eval ($f)(x::ModelType, arg::Integer) = mexpr(:call, ($f), _expr(x), arg)
+    @eval ($f)(x::ModelType, arg::Integer) = mexpr(:call, $(Expr(:quote, f)), _expr(x), arg)
 
-    @eval ($f)(x::ModelType, args...) = mexpr(:call, ($f), _expr(x), map(_expr, args)...)
+    @eval ($f)(x::ModelType, args...) = mexpr(:call, $(Expr(:quote, f)), _expr(x), map(_expr, args)...)
 end
 
 # Non-Base functions:
 for f in [:der, :pre]
-    @eval ($f)(x::ModelType, args...) = mexpr(:call, ($f), _expr(x), args...)
+    @eval ($f)(x::ModelType, args...) = mexpr(:call, $(Expr(:quote, f)), _expr(x), args...)
 end
 
 

src/sim.jl

@@ -17,12 +17,12 @@ arguments are supported (use one or the other after the first
 argument).
 
 ```julia
-sim(m::Model, tstop::Float64=1.0, Nsteps::Int=500, reltol::Float64=1e-4, abstol::Float64=1e-4)
-sim(m::Model; tstop::Float64=1.0, Nsteps::Int=500, reltol::Float64=1e-4, abstol::Float64=1e-4)
-sim(m::Sim, tstop::Float64=1.0, Nsteps::Int=500, reltol::Float64=1e-4, abstol::Float64=1e-4)
-sim(m::Sim; tstop::Float64=1.0, Nsteps::Int=500, reltol::Float64=1e-4, abstol::Float64=1e-4)
-sim(m::SimState, tstop::Float64=1.0, Nsteps::Int=500, reltol::Float64=1e-4, abstol::Float64=1e-4)
-sim(m::SimState; tstop::Float64=1.0, Nsteps::Int=500, reltol::Float64=1e-4, abstol::Float64=1e-4)
+sim(m::Model, tstop::Float64=1.0, Nsteps::Int=500, reltol::Float64=1e-4, abstol::Float64=1e-4, init::Symbol=:Ya_Ydp)
+sim(m::Model; tstop::Float64=1.0, Nsteps::Int=500, reltol::Float64=1e-4, abstol::Float64=1e-4, init::Symbol=:Ya_Ydp)
+sim(m::Sim, tstop::Float64=1.0, Nsteps::Int=500, reltol::Float64=1e-4, abstol::Float64=1e-4, init::Symbol=:Ya_Ydp)
+sim(m::Sim; tstop::Float64=1.0, Nsteps::Int=500, reltol::Float64=1e-4, abstol::Float64=1e-4, init::Symbol=:Ya_Ydp)
+sim(m::SimState, tstop::Float64=1.0, Nsteps::Int=500, reltol::Float64=1e-4, abstol::Float64=1e-4, init::Symbol=:Ya_Ydp)
+sim(m::SimState; tstop::Float64=1.0, Nsteps::Int=500, reltol::Float64=1e-4, abstol::Float64=1e-4, init::Symbol=:Ya_Ydp)
 ```
 
 ### Arguments
@@ -34,6 +34,10 @@ sim(m::SimState; tstop::Float64=1.0, Nsteps::Int=500, reltol::Float64=1e-4, abst
 * `Nsteps::Int` : the number of simulation steps, default = 500
 * `reltol::Float64` : the relative tolerance, default = 1e-4
 * `abstol::Float64` : the absolute tolerance, default = 1e-4
+* `init` : initialization of the model; options include:
+  * `:none` : no initialization
+  * `:Ya_Ydp` :  given `Y_d`, calculate `Y_a` and `Y'_d` (the default)
+  * `:Y` :  given `Y'`, calculate `Y`
 
 ### Returns
 

src/simcreation.jl

@@ -43,6 +43,8 @@ type Sim
     eq::EquationSet           # the input
     F::SimFunctions
     id::Array{Int, 1}         # indicates whether a variable is algebraic or differential
+    yfixed::Array{Bool, 1}    # indicates whether a variable is fixed
+    ypfixed::Array{Bool, 1}   # indicates whether a derivative is fixed
     outputs::Array{ASCIIString, 1} # output labels
     unknown_idx_map::Dict     # symbol => index into y (or yp)
     discrete_inputs::Set      # Discrete variables
@@ -66,8 +68,10 @@ simulation, including a Sim. Normally created with
 """ ->
 type SimState
     t::Array{Float64, 1}      # time
-    y0::Array{Float64, 1}     # state vector
-    yp0::Array{Float64, 1}    # derivatives vector
+    y0::Array{Float64, 1}     # initial state vector
+    yp0::Array{Float64, 1}    # initial derivatives vector
+    y::Array{Float64, 1}      # state vector
+    yp::Array{Float64, 1}     # derivatives vector
     structural_change::Bool
     history::SimStateHistory
     sm::Sim # reference to a Sim
@@ -103,6 +107,8 @@ function create_sim(eq::EquationSet)
 
     sm.outputs = fill_from_map("", N_unknowns, sm.y_map, x -> x.label)
     sm.id = fill_from_map(-1, N_unknowns, sm.yp_map, x -> 1)
+    sm.yfixed = fill_from_map(true, N_unknowns, sm.y_map, x -> x.fixed)
+    sm.ypfixed = fill_from_map(true, N_unknowns, sm.yp_map, x -> x.fixed)
     sm.abstol = 1e-4
     sm.reltol = 1e-4
 
@@ -135,8 +141,10 @@ function create_simstate (sm::Sim)
     N_unknowns = sm.varnum - 1
 
     t = [0.0]
-    y0 = fill_from_map(0.0, N_unknowns, sm.y_map, x -> to_real(x.value))
-    yp0 = fill_from_map(0.0, N_unknowns, sm.yp_map, x -> to_real(x.value))
+    y = fill_from_map(0.0, N_unknowns, sm.y_map, x -> to_real(x.value))
+    yp = fill_from_map(0.0, N_unknowns, sm.yp_map, x -> to_real(x.value))
+    y0 = copy(y)
+    yp0 = copy(yp)
     structural_change = false
     history = SimStateHistory (Dict(),Dict())
     for (k,v) in sm.y_map
@@ -145,7 +153,7 @@ function create_simstate (sm::Sim)
             history.x[k] = Any[]
         end
     end
-    ss = SimState (t,y0,yp0,structural_change,history,sm)
+    ss = SimState(t,y0,yp0,y,yp,structural_change,history,sm)
 
     ss
 end