Skip to content

Instantly share code, notes, and snippets.

@josejuan

josejuan/planificar-rutas.fs

Created March 26, 2014 22:46
Show Gist options
  • Save josejuan/9795379 to your computer and use it in GitHub Desktop.
Save josejuan/9795379 to your computer and use it in GitHub Desktop.
Download ZIP
Raw
planificar-rutas.fs
open Optimization
open System.Linq
type CheckPoint = {
CPId : string // checkpoint Id
Lat : double // latitude
Lon : double // longitude
}
type InterestPoint = {
IPId : string // interest point Id
CPId : string // checkpoint Id
PTId : string // (interest)point type Id
Prices : double list // prices per day
}
type UserPointTypePreferences = {
PTId : string // (interest)point type Id
MinPerCheckpoint : double // minimum interestpoint (of that type) on each checkpoints
MaxPerCheckpoint : double
MinTotal : double // minimum interestpoint (of that type) on all checkpoints
MaxTotal : double
}
type UserPreferences = {
MinDays : double // minimum total days
MaxDays : double
MinDayDistance : double // minimum checkpoint to next checkpoint distance
MaxDayDistance : double
MinTotalDistance : double // minimum travel total distance
MaxTotalDistance : double
MinTotalPrice : double // minimum travel total price
MaxTotalPrice : double
MinDaysPerCheckPoint : double // minimum days on each checkpoint
MaxDaysPerCheckPoint : double
DistanceImportance : double // how important is maximize/minimize total distance
PriceImportance : double
IPNumberImportance : double
PTPrefs : UserPointTypePreferences list
}
let checkPoints =
[ {CPId = "Ciudad A"; Lat = 43.3618742; Lon = -8.4126837}
{CPId = "Ciudad B"; Lat = 42.8802049; Lon = -8.5447697}
{CPId = "Ciudad C"; Lat = 42.2260892; Lon = -8.7254259}
{CPId = "Ciudad D"; Lat = 42.3383626; Lon = -7.8636856}
{CPId = "Ciudad E"; Lat = 43.0122837; Lon = -7.5565001}
{CPId = "Ciudad F"; Lat = 43.5705136; Lon = -7.2597601}
]
let distance (a : CheckPoint) (b : CheckPoint) : double =
sqrt ((a.Lat - b.Lat)**2. + (a.Lon - b.Lon)**2.)
let interestPoints =
[ {IPId = "Playa A" ; CPId = "Ciudad A"; PTId = "Costa" ; Prices = [ 1.; 1.; 1.; 1.; 1.; 1.; 1.; 1.; 1.; 1.]}
{IPId = "Playa B" ; CPId = "Ciudad B"; PTId = "Costa" ; Prices = [ 1.; 1.; 1.; 1.; 1.; 1.; 1.; 1.; 1.; 1.]}
{IPId = "Playa C" ; CPId = "Ciudad C"; PTId = "Costa" ; Prices = [ 1.; 1.; 1.; 1.; 1.; 1.; 1.; 1.; 1.; 1.]}
{IPId = "Playa D" ; CPId = "Ciudad D"; PTId = "Costa" ; Prices = [ 1.; 1.; 1.; 1.; 1.; 1.; 1.; 1.; 1.; 1.]}
{IPId = "Playa E" ; CPId = "Ciudad E"; PTId = "Costa" ; Prices = [ 1.; 1.; 1.; 1.; 1.; 1.; 1.; 1.; 1.; 1.]}
{IPId = "Playa F" ; CPId = "Ciudad F"; PTId = "Costa" ; Prices = [ 1.; 1.; 1.; 1.; 1.; 1.; 1.; 1.; 1.; 1.]}
{IPId = "Cultura A 1"; CPId = "Ciudad A"; PTId = "Cultura" ; Prices = [ 2.; 2.; 2.; 0.; 2.; 2.; 2.; 2.; 2.; 2.]}
{IPId = "Cultura A 2"; CPId = "Ciudad A"; PTId = "Cultura" ; Prices = [ 4.; 4.; 4.; 4.; 4.; 4.; 0.; 0.; 0.; 0.]}
{IPId = "Cultura B 1"; CPId = "Ciudad B"; PTId = "Cultura" ; Prices = [ 2.; 2.; 2.; 0.; 3.; 3.; 3.; 2.; 2.; 2.]}
{IPId = "Cultura C 1"; CPId = "Ciudad C"; PTId = "Cultura" ; Prices = [ 2.; 3.; 3.; 3.; 0.; 2.; 2.; 2.; 2.; 2.]}
{IPId = "Cultura D 1"; CPId = "Ciudad D"; PTId = "Cultura" ; Prices = [ 2.; 2.; 2.; 1.; 0.; 2.; 5.; 5.; 2.; 2.]}
{IPId = "Cultura D 2"; CPId = "Ciudad D"; PTId = "Cultura" ; Prices = [ 0.; 0.; 0.; 0.; 9.; 9.; 9.; 0.; 0.; 0.]}
{IPId = "Cultura E 1"; CPId = "Ciudad E"; PTId = "Cultura" ; Prices = [ 2.; 2.; 2.; 0.; 3.; 3.; 3.; 2.; 2.; 2.]}
{IPId = "Cultura F 1"; CPId = "Ciudad F"; PTId = "Cultura" ; Prices = [ 2.; 2.; 2.; 1.; 0.; 2.; 2.; 2.; 2.; 2.]}
{IPId = "Popular A 1"; CPId = "Ciudad A"; PTId = "Popular" ; Prices = [ 4.; 3.; 3.; 3.; 0.; 4.; 4.; 0.; 0.; 0.]}
{IPId = "Popular B 1"; CPId = "Ciudad B"; PTId = "Popular" ; Prices = [ 5.; 4.; 6.; 0.; 9.; 9.; 9.; 0.; 0.; 0.]}
{IPId = "Popular C 1"; CPId = "Ciudad C"; PTId = "Popular" ; Prices = [ 2.; 2.; 2.; 0.; 2.; 2.; 2.; 2.; 2.; 2.]}
{IPId = "Popular D 1"; CPId = "Ciudad D"; PTId = "Popular" ; Prices = [ 2.; 2.; 2.; 0.; 3.; 3.; 3.; 2.; 2.; 2.]}
{IPId = "Popular E 1"; CPId = "Ciudad E"; PTId = "Popular" ; Prices = [ 0.; 0.; 0.; 1.; 0.; 2.; 0.; 0.; 0.; 0.]}
{IPId = "Popular F 1"; CPId = "Ciudad F"; PTId = "Popular" ; Prices = [ 2.; 2.; 2.; 1.; 0.; 2.; 2.; 2.; 2.; 2.]}
{IPId = "Hotel A 1" ; CPId = "Ciudad A"; PTId = "Alojamiento"; Prices = [30.; 30.; 30.; 30.; 45.; 45.; 50.; 30.; 30.; 0.]}
{IPId = "Hotel A 2" ; CPId = "Ciudad A"; PTId = "Alojamiento"; Prices = [36.; 0.; 20.; 20.; 0.; 0.; 0.; 20.; 20.; 20.]}
{IPId = "Hotel A 3" ; CPId = "Ciudad A"; PTId = "Alojamiento"; Prices = [35.; 35.; 0.; 0.; 30.; 30.; 20.; 20.; 15.; 15.]}
{IPId = "Hotel B 1" ; CPId = "Ciudad B"; PTId = "Alojamiento"; Prices = [20.; 20.; 30.; 35.; 0.; 0.; 50.; 20.; 20.; 20.]}
{IPId = "Hotel B 2" ; CPId = "Ciudad B"; PTId = "Alojamiento"; Prices = [20.; 20.; 30.; 35.; 0.; 0.; 50.; 20.; 20.; 20.]}
{IPId = "Hotel C 1" ; CPId = "Ciudad C"; PTId = "Alojamiento"; Prices = [20.; 20.; 30.; 35.; 0.; 0.; 50.; 20.; 20.; 20.]}
{IPId = "Hotel D 1" ; CPId = "Ciudad D"; PTId = "Alojamiento"; Prices = [20.; 20.; 30.; 35.; 0.; 0.; 50.; 20.; 20.; 20.]}
{IPId = "Hotel D 2" ; CPId = "Ciudad D"; PTId = "Alojamiento"; Prices = [20.; 20.; 30.; 35.; 0.; 0.; 50.; 20.; 20.; 20.]}
{IPId = "Hotel E 1" ; CPId = "Ciudad E"; PTId = "Alojamiento"; Prices = [20.; 20.; 30.; 35.; 0.; 0.; 50.; 20.; 20.; 20.]}
{IPId = "Hotel E 2" ; CPId = "Ciudad E"; PTId = "Alojamiento"; Prices = [20.; 20.; 30.; 35.; 0.; 0.; 50.; 20.; 20.; 20.]}
{IPId = "Hotel E 3" ; CPId = "Ciudad E"; PTId = "Alojamiento"; Prices = [20.; 20.; 30.; 35.; 0.; 0.; 50.; 20.; 20.; 20.]}
{IPId = "Hotel F 1" ; CPId = "Ciudad F"; PTId = "Alojamiento"; Prices = [20.; 20.; 30.; 35.; 0.; 0.; 50.; 20.; 20.; 20.]}
]
let userPreferences = {
MinDays = 3.
MaxDays = 6.
MinDayDistance = 0.
MaxDayDistance = 0.68
MinTotalDistance = 0.
MaxTotalDistance = 1e6
MinTotalPrice = 0.
MaxTotalPrice = 999999.
MinDaysPerCheckPoint = 0. // puede fácilmente especializarse por cada checkpoint
MaxDaysPerCheckPoint = 1.
DistanceImportance = 0.1
PriceImportance = -1.
IPNumberImportance = 0.
PTPrefs =
[ {PTId = "Costa" ; MinPerCheckpoint = 0.; MaxPerCheckpoint = 9.; MinTotal = 1.; MaxTotal = 99.}
{PTId = "Cultura" ; MinPerCheckpoint = 0.; MaxPerCheckpoint = 9.; MinTotal = 1.; MaxTotal = 99.}
{PTId = "Popular" ; MinPerCheckpoint = 0.; MaxPerCheckpoint = 9.; MinTotal = 1.; MaxTotal = 99.}
{PTId = "Alojamiento"; MinPerCheckpoint = 1.; MaxPerCheckpoint = 1.; MinTotal = 0.; MaxTotal = 99.}
]
}
type CheckPointVar = {cp : CheckPoint; v : Variable}
type PriceVar = {ip : InterestPoint; v : Variable}
type EdgeVar = {a : CheckPoint; b : CheckPoint; v : Variable}
type ExpressionOpClass = ExpressionOpClass with
static member Sum(x: Variable list) = Expression.Sum x
static member Sum(x: Expression list) = Expression.Sum x
static member Sum(x: Term list) = Expression.Sum x
static member Equ(a: Variable, b: Expression) = Expression.op_Equality(a, b)
static member Equ(a: Expression, b: double) = Expression.op_Equality(a, b)
static member Leq(a: Expression, b: double) = Expression.op_LessThanOrEqual(a, b)
static member Leq(a: Expression, b: Variable) = Expression.op_LessThanOrEqual(a, b)
static member Leq(a: Variable, b: Variable) = Expression.op_LessThanOrEqual(a, b * 1.)
static member BetweenLeft(a: Expression, b: Expression) = (b, Expression.op_LessThanOrEqual(a, b))
static member BetweenLeft(a: Variable, b: Expression) = (b, Expression.op_LessThanOrEqual(a, b))
static member BetweenLeft(a: double, b: Expression) = (b, Expression.op_LessThanOrEqual(a, b))
static member BetweenRight((b, d): Expression * Constraint, c: Expression) = [d; Expression.op_LessThanOrEqual(b, c)]
static member BetweenRight((b, d): Expression * Constraint, c: double) = [d; Expression.op_LessThanOrEqual(b, c)]
// Σ xs
let inline Σ x =
let inline isum (z: ^z, x: ^x) = ((^z or ^x) : (static member Sum: ^x -> Expression) x)
isum (ExpressionOpClass, x)
// a == b
let inline (==) a b =
let inline eq (z: ^z, a: ^a, b: ^b) = ((^z or ^a or ^b) : (static member Equ: ^a * ^b -> Constraint) (a, b))
eq (ExpressionOpClass, a, b)
// a <= b
let inline (<=) a b =
let inline eq (z: ^z, a: ^a, b: ^b) = ((^z or ^a or ^b) : (static member Leq: ^a * ^b -> Constraint) (a, b))
eq (ExpressionOpClass, a, b)
// a <= b <= c
let inline (.<=) a b =
let inline bleft (z: ^z, a: ^a, b: ^b) = ((^z or ^a or ^b) : (static member BetweenLeft: ^a * ^b -> Expression * Constraint) (a, b))
bleft (ExpressionOpClass, a, b)
let inline (<=.) a b =
let inline bright (z: ^z, a: ^a, b: ^b) = ((^z or ^a or ^b) : (static member BetweenRight: ^a * ^b -> Constraint list) (a, b))
bright (ExpressionOpClass, a, b)
let solveProblem (chkpts : CheckPoint list)
(ipts : InterestPoint list)
(userp : UserPreferences ) : (int * InterestPoint list) list =
let model = new Model()
let ndays = ipts.Head.Prices.Length
let days = [0..ndays - 1]
let mkVar name args lo hi =
let v = new Variable(System.String.Format(name, List.toArray args), lo, hi, VariableType.Integer)
model.AddVariable v
v
let zero = mkVar "Z" [] 0. 0.
let binVar name args = mkVar name args 0. 1.
// Variables D1, D2, ... Dn que indican que días forman parte del trayecto (y holgura para consecutivos)
let D = [|for d in days -> binVar "D_{0}" [d]|]
let H = [|for d in [0..ndays] -> binVar "H_{0}" [d]|]
// Variables C1, C2, ... que indican para cada día, en que checkpoint se hace parada
let C = [| for d in days
-> [| for cp in chkpts
-> {cp = cp; v = binVar "C_{0}_{1}" [d; cp.CPId] } |] |]
// Variables P1, P2, ... que indican para cada precio posible de cada punto de interés de cada día si se selecciona o no
let P = [| for d in days
-> [| for ip in ipts do
if ip.Prices.[d] > 0.
then yield {ip = ip; v = binVar "P_{0}_{1}" [d; ip.IPId] } |] |]
// Variables E1, E2, ... que indican que trayecto y en que dirección se realiza cada día, desde un checkpoint al siguiente
let E = [| for d in [0..ndays - 2]
-> [| for a in [0..chkpts.Length-1] do
for b in [0..chkpts.Length-1]
-> { a = chkpts.[a]
b = chkpts.[b]
v = binVar "E_{0}_{1}_{2}" [d; chkpts.[a].CPId; chkpts.[b].CPId] } |] |]
// RESTRICCION #1
// model.AddConstraints [ for d in days -> Σ [ for cpd in C.[d] -> cpd.v ] <= 1. ]
// RESTRICCION #2
model.AddConstraints [ for d in days -> D.[d] == Σ [ for cpd in C.[d] -> cpd.v ]]
// RESTRICCION #3
model.AddConstraints (userp.MinDays .<= Σ [ for d in days -> D.[d] ] <=. userp.MaxDays)
// RESTRICCION #4
model.AddConstraints [
for cp in chkpts do
yield! let checkpointDays = [ for d in days do
for cpd in C.[d] do
if cpd.cp.CPId = cp.CPId
then yield cpd.v ] in
userp.MinDaysPerCheckPoint .<= Σ checkpointDays <=. userp.MaxDaysPerCheckPoint ]
// RESTRICCION #5
model.AddConstraints [ for d in days do
for cpd in C.[d] do
yield! let interestPoints = [ for p in P.[d] do
if p.ip.CPId = cpd.cp.CPId
then yield p.v ] in
cpd.v .<= Σ interestPoints <=. cpd.v * (double interestPoints.Length) ]
// RESTRICCION #6
model.AddConstraints [
for d in days do
for cpd in C.[d] do
for pt in userp.PTPrefs do
yield! let ip_number = Σ [ for p in P.[d] do
if p.ip.CPId = cpd.cp.CPId && p.ip.PTId = pt.PTId
then yield p.v ] in
pt.MinPerCheckpoint * cpd.v .<= ip_number <=. pt.MaxPerCheckpoint * cpd.v ]
// RESTRICCION #7
model.AddConstraints [
for pt in userp.PTPrefs do
yield! let ip_number = Σ [ for d in days do
for cpd in C.[d] do
for p in P.[d] do
if p.ip.CPId = cpd.cp.CPId && p.ip.PTId = pt.PTId
then yield p.v ] in
pt.MinTotal .<= ip_number <=. pt.MaxTotal ]
// RESTRICCION #8
let total_price = Σ [ for d in days do for p in P.[d] -> p.v * p.ip.Prices.[d] ]
model.AddConstraints ( userp.MinTotalPrice .<= total_price <=. userp.MaxTotalPrice )
// RESTRICCION #9
let uvh = [ for d in [0..ndays] do
yield ( if d = 0 then zero else D.[d - 1]
, if d = ndays then zero else D.[d]
, H.[d] ) ]
model.AddConstraints [ for (u, v, h) in uvh do yield! [ h <= u ; h <= v ; u + v - 1. <= h]]
model.AddConstraint( Σ [ for (u, v, h) in uvh -> u + v - 2. * h ] == 2. )
// Para cada par de checkpoints U y V (en cada par de días consecutivos) su arista es
// seleccionada, si los dos están seleccionados, es decir un AND:
//
// UV = U & V
//
let edges = [
for d in [0..ndays - 2] do
for uv in E.[d] do
yield let lookup w (i : CheckPoint) = C.[d + w].First(fun k -> k.cp.CPId = i.CPId).v
in (uv.v, lookup 0 uv.a, lookup 1 uv.b, distance uv.a uv.b)]
// RESTRICCION #10
model.AddConstraints [
for (ev, av, bv, dst) in edges do
yield! [ ev <= av; ev <= bv; av + bv - 1. <= ev
; userp.MinDayDistance * ev <= dst
; dst * ev <= userp.MaxDayDistance ]]
// RESTRICCIÓN #11
let total_distance = Σ [ for (ev, av, bv, dst) in edges -> ev * dst ]
model.AddConstraints (userp.MinTotalDistance .<= total_distance <=. userp.MaxTotalDistance)
// La función de optimización está parametrizada ponderando los factores a considerar
let total_ips = Σ [for p in P do for q in p -> q.v]
model.AddObjective(
total_price * userp.PriceImportance +
total_distance * userp.DistanceImportance +
total_ips * userp.IPNumberImportance, null, ObjectiveSense.Maximize)
// Exigir un determinado checkpoint inicial y otro final sería detectando fecha de arrival/departure y ligándola a su checkpoint.
let solution = (new Optimization.Solver.GLPK.GLPKSolver()).Solve(model)
let setted (v : Variable) = solution.GetVariableValue(v.Name).Value = 1.
if solution.Status <> Solver.SolutionStatus.Optimal
then []
else [for d in days do
if setted D.[d] then
yield (d, [for p in P.[d] do if setted p.v then yield p.ip ] ) ]
[<EntryPoint>]
let main argv =
let solution = solveProblem checkPoints interestPoints userPreferences
printfn "\n\n\n"
match solution with
| [] -> printfn "Not solution found!"
| xs -> do
printfn "Solution found:"
let ndays = xs.Length
let dst = [
for d in [0..ndays - 2]
-> let cp i = checkPoints.First(fun c -> c.CPId = (snd xs.[i]).Head.CPId) in
distance (cp d) (cp (d + 1))]
printfn " + Total travel price: %.2f eur" ([for (d, ips) in xs do yield! [for p in ips -> p.Prices.[d]]].Sum())
printfn " + Total travel days: %d" ndays
printfn " + Total travel distance: %f Km" (dst.Sum())
printfn " + Planning:"
[for d in [0..ndays-1] -> (d, fst xs.[d], snd xs.[d])].ToList().ForEach(
fun (i, d, ips) ->
printf " #%d Day %d (%s)" (i + 1) (d + 1) ips.Head.CPId
if i < ndays - 1
then printfn " ~ %f Km" dst.[i]
else printfn ""
ips.ToList().ForEach(
fun p ->
printfn " * %-20s: %-20s, price = %6.2f eur" p.PTId p.IPId p.Prices.[d]
)
printfn ""
)
0 // devolver un código de salida entero
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment