Merge pull request #4 from ChrisRackauckas/diffeq

Change to DifferentialEquations.jl

REQUIRE

@@ -1,3 +1,4 @@
 julia 0.5-
 MathProgBase
-ODE
+OrdinaryDiffEq
+Ipopt

docs/installation.rst

@@ -8,10 +8,10 @@ You can easily install the package by running
   Pkg.add("SwitchTimeOpt")
 
 
-This does not automatically install any nonblinear solver. To install, for example, `Ipopt <https://github.com/JuliaOpt/Ipopt.jl/>`_, just run
+This does automatically install the Ipopt nonlinear solver. To install, for example, `KNITRO <https://github.com/JuliaOpt/KNITRO.jl>`_, just run
 
 ::
     
-    Pkg.add("Ipopt")
+    Pkg.add("KNITRO")
 
 To install other solvers we refer the user to the `JuliaOpt <http://www.juliaopt.org/>`_ page.

examples/fishing.jl

@@ -15,9 +15,8 @@ plt[:rc]("font", family="serif")
 maxiter = 20
 using Ipopt
 solver = IpoptSolver(
-  print_level=0,  # Suppress output
-  max_iter = maxiter,
-  linear_solver="ma57")
+  print_level=0,
+  max_iter=maxiter)
 
 
 ### Define system parameters
@@ -139,7 +138,7 @@ end
 
 @printf("\nOptimal Switching Times for ngrid = 200\n")
 @printf("--------------------------------------\n")
-@printf("tauopt = "); show(round(tauopt[:, idx200],3)); @printf("\n")
+@printf("tauopt = "); show(round.(tauopt[:, idx200],3)); @printf("\n")
 
 
 # Generate plots for ngrid = 25

examples/linear.jl

@@ -12,9 +12,10 @@ plt[:rc]("text", usetex=true)
 plt[:rc]("font", family="serif")   # Use Serif math
 
 using Ipopt
-solver = IpoptSolver(
-          print_level = 0,
-          linear_solver="ma57")
+# solver = IpoptSolver(
+#           print_level = 0,
+#           linear_solver="ma57")
+solver = IpoptSolver(print_level = 0)
 
 
 ### Define system parameters
@@ -66,7 +67,7 @@ xsim, xpts, Jsim, t = simulate(m)
 @printf("Objective Function at the Optimum:              J = %.3f\n", Jopt)
 @printf("Simulated Objective Function at the Optimum: Jsim = %.3f\n", Jsim)
 @printf("Elapsed Time:                                time = %.2f [ms]\n", mean(soltime)*10^3)
-@printf("Optimum:                                      tau = "); show(round(tauopt,3)); @printf("\n")
+@printf("Optimum:                                      tau = "); show(round.(tauopt,3)); @printf("\n")
 
 ### Plot results
 figure()

examples/tank.jl

@@ -19,8 +19,7 @@ using Ipopt
 solver = IpoptSolver(
   print_level=0,
   tol = 1e-10,  # Increase required precision to avoid convergence before maxiter iterations
-  max_iter = maxiter,
-  linear_solver="ma57")
+  max_iter = maxiter)
 
 
 ### Define system parameters
@@ -143,7 +142,7 @@ end
 
 @printf("\nOptimal Switching Times for ngrid = 10\n")
 @printf("--------------------------------------\n")
-@printf("tauopt = "); show(round(tauopt[:, 1],2)); @printf("\n")
+@printf("tauopt = "); show(round.(tauopt[:, 1],2)); @printf("\n")
 
 
 

src/SwitchTimeOpt.jl

@@ -2,7 +2,9 @@ module SwitchTimeOpt
 
   # Import Necessary Modules
   using MathProgBase
-  using ODE
+  using DiffEqBase
+  using OrdinaryDiffEq
+  using Ipopt
 
 
   export stoproblem,

src/interface.jl

@@ -10,7 +10,7 @@ function stoproblem(
   lb::Array{Float64, 1}=emptyfvec,      # Lower Bound on Intervals
   ub::Array{Float64, 1}=emptyfvec,      # Upper Bound on Intervals
   tau0ws::Array{Float64,1}=emptyfvec,   # Warm Start tau vector
-  solver::MathProgBase.AbstractMathProgSolver=MathProgBase.defaultNLPsolver)
+  solver::MathProgBase.AbstractMathProgSolver=Ipopt.IpoptSolver())
 
   # Get Dimensions
   nx = size(A, 1)      # State Dimension
@@ -141,7 +141,7 @@ function stoproblem(
   lb::Array{Float64, 1}=emptyfvec,  # Lower Bound on Intervals
   ub::Array{Float64, 1}=emptyfvec,  # Upper Bound on Intervals
   tau0ws::Array{Float64,1}=emptyfvec, # Warm Start tau vector
-  solver::MathProgBase.AbstractMathProgSolver=MathProgBase.defaultNLPsolver)
+  solver::MathProgBase.AbstractMathProgSolver=Ipopt.IpoptSolver())
 
   # Get Dimensions
   nx = length(x0)         # State Dimension

src/simulation.jl

@@ -298,16 +298,16 @@ function simulatenlinsto(nonlin_dyn::Function, tau::Array{Float64,1}, x0::Array{
     end
 
     if tempInd2>tempInd1  # There has been a progress and the switching times are not collapsed. So we integrate.
+      prob = DiffEqBase.ODEProblem(nldyn,xprevSwitch,(t[tempInd1],t[tempInd2]))
       if tempInd2 == tempInd1 + 1  # Only one step progress. Handle as special case for the integrator
-
-        _, xmap = ode45(nldyn, xprevSwitch, [t[tempInd1]; (t[tempInd1]+t[tempInd2])/2; t[tempInd2]], points=:specified)  # Integrate
-        xmap = hcat(xmap...)
-        xtemp = [xmap[:,1] xmap[:,3]]          # Take only two points
+        sol = DiffEqBase.solve(prob,OrdinaryDiffEq.Tsit5();save_everystep=false)
+        xtemp = [sol[1] sol[end]]
+        #_, xmap = ode45(nldyn, xprevSwitch, [t[tempInd1]; (t[tempInd1]+t[tempInd2])/2; t[tempInd2]], points=:specified)  # Integrate
+        #xmap = hcat(xmap...)
+        #xtemp = [xmap[:,1] xmap[:,3]]          # Take only two points
       else
-
-        _, xmap = ode45(nldyn, xprevSwitch, t[tempInd1:tempInd2], points=:specified)
-        xtemp = hcat(xmap...)  # https://github.com/JuliaLang/ODE.jl/issues/80
-
+        sol = DiffEqBase.solve(prob,OrdinaryDiffEq.Tsit5();saveat=t[tempInd1:tempInd2])
+        xtemp = hcat((sol[i] for i in eachindex(sol))...)
       end
 
       x[:, tempInd1:tempInd2] = xtemp