added tests for benchmarking and removed @parallel because it was messing things up

Author: huckl3b3rry87

REQUIRE

@@ -8,3 +8,4 @@ KNITRO
 Interpolations
 MathProgBase
 CSV
+Images

src/Base.jl

@@ -173,17 +173,17 @@ function interpolateLagrange!(n; numPts::Int64=250, tfOptimal::Any=false)
     x_int, u_int = intervals(n, int, copy(n.r.X), n.r.U)
 
     # sample polynomial in interval at n.r.t_polyPts
-    @parallel for st in 1:n.numStates
+     for st in 1:n.numStates
       n.r.X_polyPts[st][int] = interpolate_lagrange(n.r.t_polyPts[int], t_st_int[int], x_int[:,st])'
     end
 
-    @parallel for ctr in 1:n.numControls
+     for ctr in 1:n.numControls
       n.r.U_polyPts[ctr][int] = interpolate_lagrange(n.r.t_polyPts[int], t_ctr_int[int], u_int[:,ctr])'
     end
 
     # sample polynomial in interval at n.r.t_polyPts NOTE costate is missing the last point, that is the t_st_int[int][1:end-1]
     if n.s.evalCostates && n.s.evalConstraints
-      @parallel for st in 1:n.numStates
+       for st in 1:n.numStates
         n.r.CS_polyPts[st][int] = interpolate_lagrange(n.r.t_polyPts[int], t_st_int[int][1:end-1], n.r.CS[st][int])'
       end
     end
@@ -195,20 +195,20 @@ function interpolateLagrange!(n; numPts::Int64=250, tfOptimal::Any=false)
   totalPts = length(n.r.t_pts);
 
   n.r.X_pts = Matrix{Float64}(totalPts, n.numStates)
-  @parallel for st in 1:n.numStates # states
+   for st in 1:n.numStates # states
     temp = [n.r.X_polyPts[st][int][1:end,:] for int in 1:n.Ni];
     n.r.X_pts[:,st] = [idx for tempM in temp for idx=tempM];
   end
 
   n.r.U_pts = Matrix{Float64}(totalPts, n.numControls)
-  @parallel for ctr in 1:n.numControls # controls
+   for ctr in 1:n.numControls # controls
     temp = [n.r.U_polyPts[ctr][int][1:end,:] for int in 1:n.Ni];
     n.r.U_pts[:,ctr] = [idx for tempM in temp for idx=tempM];
   end
 
   if n.s.evalCostates && n.s.evalConstraints
     n.r.CS_pts = Matrix{Float64}(totalPts, n.numStates)
-    @parallel for st in 1:n.numStates # states
+     for st in 1:n.numStates # states
       temp = [n.r.CS_polyPts[st][int][1:end,:] for int in 1:n.Ni];
       n.r.CS_pts[:,st] = [idx for tempM in temp for idx=tempM];
     end
@@ -241,12 +241,12 @@ function interpolateLinear!(n; numPts::Int64=250, tfOptimal::Any=false)
   n.r.U_pts = Matrix{Float64}(numPts, n.numControls)
 
   knots = (n.r.t_st,)
-  @parallel for st in 1:n.numStates
+  for st in 1:n.numStates
     sp_st = interpolate(knots,n.r.X[:,st],Gridded(Linear()))
     n.r.X_pts[:,st] = sp_st[n.r.t_pts]
   end
 
-  @parallel for ctr in 1:n.numControls
+   for ctr in 1:n.numControls
     sp_ctr = interpolate(knots,n.r.U[:,ctr],Gridded(Linear()))
     n.r.U_pts[:,ctr] = sp_ctr[n.r.t_pts]
   end
@@ -277,11 +277,11 @@ function plant2dfs!(n,sol)
   dfs_plant=DataFrame();
   dfs_plant[:t]=t_sample;
 
-  @parallel for st in 1:n.numStates
+  for st in 1:n.numStates
     dfs_plant[n.state.name[st]]=[sol(t)[st] for t in t_sample];
   end
 
-  @parallel for ctr in 1:n.numControls
+  for ctr in 1:n.numControls
     dfs_plant[n.control.name[ctr]]= n.r.U[ctr];
   end
 
@@ -305,10 +305,10 @@ Date Create: 2/10/2017, Last Modified: 11/10/2017 \n
 function dvs2dfs(n)
   dfs=DataFrame()
   dfs[:t]=n.r.t_st + n.mpc.t0;
-  @parallel for st in 1:n.numStates
+  for st in 1:n.numStates
     dfs[n.state.name[st]]=n.r.X[:,st];
   end
-  @parallel for ctr in 1:n.numControls
+  for ctr in 1:n.numControls
     if n.s.integrationMethod==:tm
       dfs[n.control.name[ctr]]=n.r.U[:,ctr];
     else
@@ -318,12 +318,12 @@ function dvs2dfs(n)
 
   if n.s.evalCostates && n.s.integrationMethod == :ps && n.s.evalConstraints
     CS_vector = Matrix{Float64}(n.numStatePoints, n.numStates)
-    @parallel for st in 1:n.numStates # states
+    for st in 1:n.numStates # states
       temp = [n.r.CS[st][int][1:end,:] for int in 1:n.Ni];
       CS_vector[1:end-1,st] = [idx for tempM in temp for idx=tempM];
       CS_vector[end,st] = NaN;
     end
-    @parallel for st in 1:n.numStates # states
+    for st in 1:n.numStates # states
       dfs[Symbol(n.state.name[st],:_cs)]=CS_vector[:,st];
     end
   end

src/NLOptControl.jl

@@ -10,6 +10,7 @@ using FastGaussQuadrature
 using DataFrames
 using Ranges
 using CSV
+using Images
 
 include("Base.jl")
 using .Base

test/ocp.jl

@@ -36,3 +36,100 @@ dx=Array{Expr}(2);dx[1]=:(x[j]);dx[2]=:(u[j]-1.625);
   @show n.r.dfs_opt[1][:t_solve][1]
   @test isapprox(8.9253,n.r.obj_val[1],atol=tol)
 end
+
+############################
+# Benchmarking Test
+v0 = -2
+h0 = 10
+pts = 10000
+tf_star = 2/3*v0 +4/3*(1/2*v0^2+3/2*h0)^0.5
+ts = tf_star/2+v0/3
+
+function optimal_solution(t)
+  if t < ts
+    h = -3/4*t^2 + v0*t + h0
+    v = -3/2*t + v0
+    u = 0
+  else
+    h = 3/4*t^2 + (-3*ts + v0)*t + 3/2*ts^2 + h0
+    v = 3/2*t + (-3*ts + v0)
+    u = 3
+  end
+  return h, v, u
+end
+
+t_opt = linspace(0,tf_star,pts)
+h_opt = zeros(pts); v_opt = zeros(pts); u_opt = zeros(pts);
+for i in 1:pts
+  h, v, u = optimal_solution(t_opt[i])
+  h_opt[i] = h; v_opt[i] = v; u_opt[i] = u;
+end
+
+opt_runs = 2
+col_pts = 10:2:12
+Nck_vec = [[col_pts[i]] for i in 1:length(col_pts)]
+opt_num = length(Nck_vec)
+
+@testset "MoonLander test to make sure that a benchmark test will be able to be executed with (:integrationScheme=>$(integrationConfig)) " for integrationConfig in integrationConfigs
+
+  # final average results
+  t_opt_ave = NaN*zeros(opt_num)
+  h_error_ave = NaN*zeros(opt_num)
+  v_error_ave = NaN*zeros(opt_num)
+  u_error_ave = NaN*zeros(opt_num)
+  max_error_ave = NaN*zeros(opt_num)
+
+  for num in 1:length(col_pts)
+    n=define(numStates=2,numControls=1,X0=[h0,v0],XF=[0.,0.],XL=[-20,-20],XU=[20,20],CL=[0.],CU=[3.]);
+    n.s.tf_max = 1000.0
+    n.s.numInterpPts = pts  # must be the same as above to calculate error
+    n.s.tfOptimal = tf_star # must be the same as above to calculate error
+    states!(n,[:h,:v];descriptions=["h(t)","v(t)"]);
+    controls!(n,[:T];descriptions=["T(t)"]);
+    dx=[:(v[j]),:(T[j]-1.5)]
+    dynamics!(n,dx)
+
+    if (integrationConfig == :lgrImplicit || integrationConfig == :lgrExplicit)
+      Nck = Nck_vec[num]
+      configure!(n;(:integrationScheme=>integrationConfig),(:Nck=>Nck),(:finalTimeDV=>true));
+    else
+      configure!(n;(:integrationScheme=>integrationConfig),(:N=>col_pts[num]),(:finalTimeDV=>true));
+    end
+    x1=n.r.x[:,1];x2=n.r.x[:,2];
+    obj=integrate!(n,:(T[j]));
+    @NLobjective(n.mdl, Min, obj);
+    setvalue(n.tf, 1.5)
+    for i in 1:length(x1); setvalue(x1[i], 0.0); setvalue(x2[i], 0.0);  end
+    # cache functions; inital optimization
+    optimize!(n);
+
+    # RESET temp variables for averaging results
+    t_solve = NaN*zeros(opt_runs)
+    h_error = NaN*zeros(opt_runs)
+    v_error = NaN*zeros(opt_runs)
+    u_error = NaN*zeros(opt_runs)
+    max_error = NaN*zeros(opt_runs)
+
+    for j in 1:opt_runs
+      optimize!(n);
+      if n.r.status == :Optimal
+        t_solve[j] = n.r.t_solve
+        h_error[j] = maximum(abs.(n.r.X_pts[:,1] - h_opt))
+        v_error[j] = maximum(abs.(n.r.X_pts[:,2] - v_opt))
+        u_error[j] = maximum(abs.(n.r.U_pts[:,1] - u_opt))
+        max_error[j] = maximum([h_error[j], v_error[j]])
+      end
+    end
+
+    t_opt_ave[num] = Images.meanfinite(t_solve,1)[1]
+    h_error_ave[num] = Images.meanfinite(h_error,1)[1]
+    v_error_ave[num] = Images.meanfinite(v_error,1)[1]
+    u_error_ave[num] = Images.meanfinite(u_error,1)[1]
+    max_error_ave[num] = Images.meanfinite(max_error,1)[1]
+  end
+  @show maximum(max_error_ave)
+  @test isapprox(0, maximum(max_error_ave),atol=big_tol)
+end
+
+# Benchmarking Test
+############################

test/runtests.jl

@@ -2,8 +2,9 @@ using NLOptControl
 using Base.Test
 
 # constants
-const tol=5e-2;
-const integrationConfigs=[:lgrImplicit,:lgrExplicit,:trapezoidal,:bkwEuler]
+const tol = 5e-2
+const big_tol = 0.5 # can reduce if the number of points are increased
+const integrationConfigs = [:lgrExplicit,:lgrImplicit,:trapezoidal,:bkwEuler]
 
 include("ocp.jl")