using PyCall
PyCall.pyversion
@pyimport runGP
include("MIGP.jl")
This notebook extends the MIGP case to a more interesting problem, the design of an electric vertical takeoff and landing aircraft. The modeling of the vehicle is in the evtol.py file - it uses first-principles models to size an eVTOL aircraft for a 60 nmi mission carrying 2 passengers. The key models are:
-Beam bending models for wing spar sizing
-Parabolic drag polar for induced drag
-Airfoil data fit for wing profile drag
-Fuselage/tail profile drag based on wetted area
-Constant specific energy batteries (300 Wh/kg)
-Actuator disk theory propellor models
-Motor models based on equivalent circuit for DC brushless motor
These models are similar to those described in this reference:
Michael J. Burton, Christopher Courtin, et.al. "Solar Aircraft Design Trade Studies Using Geometric Programming", 2018 Multidisciplinary Analysis and Optimization Conference, AIAA AVIATION Forum, (AIAA 2018-3740)
-Vehicle size is constrained to have a 50 ft wingspan, and the maximum propeller disk area must fit within a 50x50 box.
This vehicle has two sets of motors; one for vertical lift, and one for forward thrust. We'd like to leave the number of motors in each group as a free variable, and constrain it to be an integer using the branch-and-bound method. The lift motors must be supported by an integer number of support booms. Therefore, there are three integer variables in the model: Number lifting motors, number of thrust motors, and number of support booms.
function branch(cur_val,g,best_val, cleaned_vars,direction, pkg_name, mod_name)
(above, below) = nearest_ints(cur_val)
Above = g(pkg_name, mod_name,cleaned_vars, [best_val...,above])
Below = g(pkg_name, mod_name,cleaned_vars, [best_val...,below])
(val, cost, I_vr, I_vl, direct) = pick_next((Above...,above), (Below...,below))
return (val, cost, I_vr, I_vl, direct)
end
pkg_name = "evtol"
mod_name = "Mission"
t = @elapsed @time out = solve_MIGP(runGP.run_GP, pkg_name, mod_name)
println(string(out[4])*" calls to solver ("*string(t/out[4])*" sec/call)")
thrust_i = findfirst(x->occursin("thrust",x), out[1])
lift_i = findfirst(x->occursin("lift",x), out[1])
boom_i = findfirst(x->occursin("boom",x), out[1])
println("Vehicle has "*string(out[2][lift_i][1])*" lift motors and "*string(out[2][thrust_i][1])*" thrust motors")
println("There are "*string(out[2][boom_i][1])* " support booms")
println("Total vehicle weight: "*string(out[3])*" lbf")
@time out = runGP.run_GP(pkg_name,mod_name)
println("Relaxed cost:"*string(out[3])*" lbf")
These results are interesting. It shows that the cost in this case to make these variables integer is very close to the optimal relaxed solution. It also suggests that you'd like to use a relatively small number of large lifting motors (to keep the drag down from the booms, while still getting a low disk loading), and that you'd like to use a distributed propulsion scheme for the forward thrust motors, possible to get soem weight savings via lighter, lower-torque motors.
using Distributed
addprocs(2)
@everywhere using PyCall
@everywhere @pyimport runGP
function branch(cur_val,g,best_val, cleaned_vars, direction, pkg_name, mod_name)
(above, below) = nearest_ints(cur_val)
input_list = [(pkg_name, mod_name,cleaned_vars, [best_val...,above], [direction...,">="]),
(pkg_name, mod_name,cleaned_vars, [best_val...,below], [direction...,"<="])]
results = @distributed vcat for inp in input_list
g(inp...)
end
(val, cost, I_vr, I_vl, dir) = pick_next((results[1]...,above), (results[2]...,below))
return (val, cost, I_vr, I_vl, dir)
end
t = @elapsed @time out = solve_MIGP(runGP.run_GP,pkg_name, mod_name)
println(string(out[4])*" calls to solver ("*string(t/out[4])*" sec/call)")
Interestingly there is some uncertainty in how many class to the solver are required to solve the problem. This is likely due to the very low local sensivity to the decision variables at the optimum point (1 optimally designed thrust motor is very similar to two optimally designed thrust motors from an overall system point of view, but 4 is clearly worse.) But looking at the total time normalized by the number of calls to the solver (or average time per call), we do see about a 25% speedup from parallelization.
@everywhere using Dagger
function branch(cur_val,g,best_val, cleaned_vars, direction, pkg_name, mod_name)
(above, below) = nearest_ints(cur_val)
results = Dict()
input_list = [(pkg_name, mod_name,cleaned_vars, [best_val...,above], [direction...,">="]),
(pkg_name, mod_name,cleaned_vars, [best_val...,below], [direction...,"<="])]
results = collect(delayed(vcat)(delayed(g)(input_list[1]...), delayed(g)(input_list[2]...)))
(val, cost, I_vr, I_vl, dir) = pick_next((results[1]...,above), (results[2]...,below))
return (val, cost, I_vr, I_vl, dir)
end
t = @elapsed @time solve_MIGP(runGP.run_GP,pkg_name, mod_name)
println(string(out[4])*" calls to solver ("*string(t/out[4])*" sec/call)")
Alternative parallel implementations give roughly the same performance.
The thing that takes the most time with this algorithm is checking the solution is optimal after the initial assignment of all variables. When the first integer variable assignment is made, all subsequent decision variables are still free to be non-integer. However, once all subsequent variables are given an integer value the optimal choice for the first value may be different; so this algoirthm goes back and checks the optimal integer assignment variables under the constraint that all other variables are constrained. This is done iteratively, by relaxing each variable succesively and picking the new integer optimum until no changes need to be made. If this process isn't completed, the solution doesn't change much (the total cost almost not at all) and the algorithm speeds up considerably.
t = @elapsed @time out = solve_MIGP(runGP.run_GP,pkg_name, mod_name)
println(string(out[4])*" calls to solver ("*string(t/out[4])*" sec/call)")
thrust_ind = findall(x->occursin("thrust",x), out[1])
lift_ind = findall(x->occursin("lift",x), out[1])
println("Results:"*string(out[2][lift_ind][1])*" lift motors, and "*string(out[2][thrust_ind][1])*" thrust motors. Total weight: "*string(out[3])*" lbs")
t = @elapsed @time out = solve_MIGP(runGP.run_GP,pkg_name, mod_name, false)
println(string(out[4])*" calls to solver ("*string(t/out[4])*" sec/call)")
thrust_ind = findall(x->occursin("thrust",x), out[1])
lift_ind = findall(x->occursin("lift",x), out[1])
println("Results:"*string(out[2][lift_ind][1])*" lift motors, and "*string(out[2][thrust_ind][1])*" thrust motors. Total weight: "*string(out[3])*" lbs")
As can be seen, the integer variables are slightly different but the total cost function is almost identical. Practically we may be able to dispense with this process, or implement limits on the required change in objective function.
One other way we might be able to speed this up is by using inequality constraints instead of equality constraints when we branch the tree. For example, given a variable $x_1$ which has a value of .5 in the relaxed solution. The above implementation solves two branches off the relaxed problem, which add the equality constraints $x_1 = 0$ and $x_1 = 1$. However, it may be more robust to say $x_1 <= 0$ and $x_1 >= 1$, both to get faster convergence and to allow x_1 more freedom to change based on subsequent choices. The con_GP function implements this change.
include("MIGP.jl")
t = @elapsed @time out = solve_MIGP(runGP.con_GP, pkg_name, mod_name)
println(string(out[4])*" calls to solver, ("*string(t/out[4])*" sec/call)")
thrust_ind = findall(x->occursin("thrust",x), out[1])
lift_ind = findall(x->occursin("lift",x), out[1])
println("Results:"*string(out[2][lift_ind][1])*" lift motors, and "*string(out[2][thrust_ind][1])*" thrust motors. Total weight: "*string(out[3])*" lbs")
In practice this uses more solver calls to get the same result as our normal solution checking process, so this doesn't appear to have much value, at least in this example.
Future work: There are two main areas of continued development on this project that I plan to continue over the next several months. The first is in increasing the fidelity, accuracy, and complexity of the underlying aircraft model. A third type of motor (rotating to provide lift and thrust) will be added; but currently that seems to require a more complex, and slower, formulation of the problem. It's also purely in Python so I focused less on that for this project. The other area of improvement is working on how to handle binary decision variables, which may or may not have logical implications for other decision variables. For example, we would like to say that our vehicle must have an integer number of wings or rotors but must have one or the other. This can be done using the continuous integer formulation developed here if combinations of the two are admissable (the vehicle can have wings and rotors), but if they are exlusionary constraints (wings or rotors but not both) then there needs to be an additional level of logic. It's not yet clear if that's useful capability in solving this problem, or how precisely that would be implemented.