Add entropic risk measure (#347)

* Add the entropic risk measure

* Add docs

docs/src/guides/add_a_risk_measure.md

@@ -304,3 +304,32 @@ round.(risk_adjusted_probability, digits = 1)
  0.8
  0.0
 ```
+
+### Entropic
+
+```@docs
+SDDP.Entropic
+```
+
+```jldoctest intermediate_risk
+risk_measure = SDDP.Entropic(0.1)
+
+SDDP.adjust_probability(
+    risk_measure,
+    risk_adjusted_probability,
+    nominal_probability,
+    noise_supports,
+    cost_realizations,
+    is_minimization
+)
+
+round.(risk_adjusted_probability, digits = 4)
+
+# output
+
+4-element Array{Float64,1}:
+ 0.11
+ 0.1991
+ 0.3648
+ 0.326
+```

src/plugins/risk_measures.jl

@@ -541,3 +541,46 @@ function adjust_probability(
     copyto!(risk_adjusted_probability, JuMP.value.(p))
     return 0.0
 end
+
+# ================================= Entropic ============================== #
+
+"""
+    Entropic(γ::Float64)
+
+The entropic risk measure as described by Dowson, Morton, and Pagnoncelli
+(2020). Multistage stochastic programs with the entropic risk measure.
+http://www.optimization-online.org/DB_HTML/2020/08/7984.html
+"""
+mutable struct Entropic <: AbstractRiskMeasure
+    γ::Float64
+end
+
+function Base.show(io::IO, ent::Entropic)
+    return print(io, "Entropic risk measure with γ = $(ent.γ)")
+end
+
+function adjust_probability(
+    measure::Entropic,
+    Q::Vector{Float64},
+    p::Vector{Float64},
+    ::Vector,
+    X::Vector{Float64},
+    is_min::Bool,
+)
+    if measure.γ == 0.0  # Special case of entropic: if γ = 0, F[X] = E[X].
+        Q .= p
+        return 0.0
+    end
+    # Handle maximization problems by negating γ. Usually we would negate X, but
+    # with convex RM's, we also have to negate the α(q) term, so it's easier to
+    # just negate γ.
+    γ = (is_min ? 1.0 : -1.0) * measure.γ
+    # Use `BigFloat` to avoid overflow that occurs when calculating `exp(x)`.
+    y = p .* exp.(big.(γ .* X))
+    Q .= y / sum(y)
+    α = sum(
+        # To avoid numerical issues with the log, skip elements that are `≈ 0`.
+        qi * log(qi / pi) for (pi, qi) in zip(p, Q) if pi > 1e-14  && qi > 1e-14
+    )
+    return -α / γ
+end

test/plugins/risk_measures.jl

@@ -211,7 +211,7 @@ end
         )
         @test risk_adjusted_probability ≈ [0.1, 0.2, 0.3, 0.2, 0.2] atol = 1e-6
     end
-   
+
     @testset "Non-uniform Worst Case" begin
         risk_adjusted_probability = Vector{Float64}(undef, 5)
         SDDP.adjust_probability(
@@ -224,7 +224,7 @@ end
         )
         @test risk_adjusted_probability ≈ [0.0, 1.0, 0.0, 0.0, 0.0] atol = 1e-6
     end
-   
+
     @testset "Non-uniform Worst Case Max" begin
         risk_adjusted_probability = Vector{Float64}(undef, 5)
         SDDP.adjust_probability(
@@ -237,7 +237,7 @@ end
         )
         @test risk_adjusted_probability ≈ [0.0, 0.0, 0.0, 0.0, 1.0] atol = 1e-6
     end
-   
+
     @testset "Non-uniform" begin
         risk_adjusted_probability = Vector{Float64}(undef, 5)
         SDDP.adjust_probability(
@@ -410,3 +410,40 @@ end
         @test risk_adjusted_probability ≈ [0.0, 1 / 6, 5 / 6, 0.0]
     end
 end
+
+@testset "Entropic" begin
+    @test sprint(show, SDDP.Entropic(0.1)) == "Entropic risk measure with γ = 0.1"
+    @testset "Minimization" begin
+        # Test that increasing values of θ lead to larger values for F[X].
+        X = [1.0, 2.0, 3.0]
+        P = [0.5, 0.5, 0.0]
+        Q = [NaN, NaN, NaN]
+        last, last_q2 = -Inf, 0.0
+        for i = -4:4
+            θ = 10.0^i
+            α = SDDP.adjust_probability(SDDP.Entropic(θ), Q, P, [], X, true)
+            current = Q' * X + α
+            @test current >= last
+            @test Q[2] >= last_q2
+            last, last_q2 = current, Q[2]
+        end
+    end
+    @testset "Maximization" begin
+        # Test that increasing values of θ lead to smaller values for F[X].
+        X = [1.0, 2.0, 3.0]
+        P = [0.5, 0.5, 0.0]
+        Q = [NaN, NaN, NaN]
+        last, last_q1 = Inf, 0.0
+        for i = -4:4
+            θ = 10.0^i
+            α = SDDP.adjust_probability(SDDP.Entropic(θ), Q, P, [], X, false)
+            current = Q' * X + α
+            if Q[1] < 1 - 1e-6
+                # If Q[1] ≈ 1, this test can fail due to numerical error.
+                @test current <= last
+            end
+            @test Q[1] >= last_q1
+            last, last_q1 = current, Q[1]
+        end
+    end
+end