PSORLab / McCormick.jl

A forward McCormick operator library
MIT License
15 stars 7 forks source link

Issue with the subtangent of differentiable McCormick relaxations #44

Closed ca0h closed 3 years ago

ca0h commented 3 years ago

Observed an incorrect subtangent of differentiable McCormick relaxations:

MWE

Code to reproduce this figure:

using McCormick, DataFrames, IntervalArithmetic

f(x,p) = (x-p)^2

# create a data frame to store output data
df = DataFrame(x = Float64[], y = Float64[], cv1 = Float64[],
               cc1 = Float64[],  l1 = Float64[],
               u1 = Float64[], scv1 = Float64[], scc1 = Float64[])

# generate relaxations
MCOption = Diff # NS
p = 0.1
p_I = Interval(0.0,2.0)
p_mc = MC{2,MCOption}(p,p_I,2)

X = Interval(0.0,3.0)
x_range = range(X.lo,stop=X.hi,length=100)
for x in x_range
    y = f(x,p)
    x_mc = MC{2,MCOption}(x,X,1)
    f_mc = f(x_mc, p_mc)
    save_tuple = (x, y, f_mc.cv, f_mc.cc,
                   lo(f_mc.Intv), hi(f_mc.Intv),
                   f_mc.cv_grad[1], f_mc.cc_grad[1])
    push!(df, save_tuple)
end

# generate subtangent lines at posi
df2 = DataFrame(x = Float64[], tan_cv = Float64[], tan_cc = Float64[])
posi = 40
x_range2 = range(x_range[posi-10], stop=x_range[posi+10], length=20)
for x in x_range2
    tanLo = df.scv1[posi]*(x - df.x[posi]) + df.cv1[posi]
    tanHi = df.scc1[posi]*(x - df.x[posi]) + df.cc1[posi]
    push!(df2, (x, tanLo, tanHi))
end

using Plots
pyplot()
data_mc = [df.y,df.cv1,df.cc1]
plt = Plots.plot(x_range, data_mc, title="McCormick Relaxation")
plt = Plots.plot!(x_range2, df2.tan_cv, linewidth=3, label = "subtangents")
plt = Plots.plot!(x_range2, df2.tan_cc, linewidth=3, label = "supertangents")
display(plt)

Relevant environment information: versioninfo()

Julia Version 1.5.3
Commit 788b2c77c1 (2020-11-09 13:37 UTC)
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: AMD Ryzen 5 2600X Six-Core Processor
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-9.0.1 (ORCJIT, znver1)
Environment:
  JULIA_EDITOR = "C:\Users\user\AppData\Local\atom\app-1.54.0\atom.exe"  -a
  JULIA_NUM_THREADS = 12

Pkg.status()

[c52e3926] Atom v0.12.30
  [6e4b80f9] BenchmarkTools v0.5.0
  [a076750e] CPLEX v0.7.6
  [336ed68f] CSV v0.8.2
  [5ae59095] Colors v0.12.6
  [a93c6f00] DataFrames v0.22.2
  [864edb3b] DataStructures v0.17.20
  [41bf760c] DiffEqSensitivity v6.35.0
  [0c46a032] DifferentialEquations v6.16.0
  [31c24e10] Distributions v0.23.8
  [bb8be931] EAGO v0.5.2
  [f6369f11] ForwardDiff v0.10.14
  [7073ff75] IJulia v1.23.1
  [d1acc4aa] IntervalArithmetic v0.17.6
  [b6b21f68] Ipopt v0.6.5
  [4076af6c] JuMP v0.21.4
  [e5e0dc1b] Juno v0.8.4
  [53c679d3] McCormick v0.10.3
  [1dea7af3] OrdinaryDiffEq v5.42.3
  [91a5bcdd] Plots v1.10.1
  [d330b81b] PyPlot v2.9.0
  [90137ffa] StaticArrays v0.12.5
mewilhel commented 3 years ago

Thanks! Good catch! I apparently had the following in line 32 of power.jl:

((xL < 0.0) && (0.0 <= xU) && (0.0 <= x)) && (3.0*x^2)/xU

when there should have been a return statement:

((xL < 0.0) && (0.0 <= xU) && (0.0 <= x)) && return (3.0*x^2)/xU

I've fixed this in the master and I'll tag a new minor release some time tomorrow.

mewilhel commented 3 years ago

Release for package with bug fix now underway. I'm closing this so the TagBot history shows this is resolved in that release.