Non-Differentiable Constraints

Specific Domains Optimization (Mathematical)
1

I am trying to write an optimization routine for insurance book. The target is to be able to increase the return without affecting the tail value at risk by modifying the shares for each contract.
This is an example of what I had so far.

#simulated losses for each contract, for every year of simulation
losses_matrix = collect([
        93.043 51.826 1.063
        84.585 52.823 0.938
        111.453 56.477 1.000
        99.525 49.805 0.938
        95.819 50.287 1.438
        114.708 51.521 1.700
        111.515 51.531 2.540
        113.211 48.664 2.390
        104.942 55.744 3.120
        99.827 47.916 2.980
    ]')

shares = [0.5, 0.3, 0.2] #initial shares
minimum_shares = [0.1, 0.1, 0.1] #boundaries for the optimizer
maximum_shares = [0.9, 0.5, 0.4]
premiums = [160, 80.0, 10.0]  #contracts premium
rp = 2 # var point for constraint
initial_return = portfolio_return(shares)
portfolio_return(x) = sum(premiums .* x .- mean_losses(x))
objective_prem_vs_initial(x) = sum(premiums .* x .- mean_losses(x)) / initial_return

function tvar(x)
    portfolio_losses = sum(losses_matrix .* x, dims=1)  # Sum across contracts (rows)
    sorted_portfolio_losses = sort(portfolio_losses[:], rev=true)  # Flatten & sort descending
    tvar_ret = mean(sorted_portfolio_losses[1:rp])  # Take top rp values
    return tvar_ret
end

tvar_lower_bound = tvar(shares) * .8
tvar_upper_bound = tvar(shares) * 1.2

model = Model(Ipopt.Optimizer)
set_silent(model)
@variable(model,  minimum_shares[i] <= x[i in 1:num_contracts] <= maximum_shares[i])
set_start_value.(x,  shares)
#@constraint(model, tvar(x) <= tvar_upper_bound)
#@constraint(model, tvar(x) >= tvar_lower_bound)
@objective(model, Max, objective_prem_vs_initial(x))
optimize!(model)
is_success = is_solved_and_feasible(model)

The problem starts when I try to add the constraint above, which I understand it cannot work as the sort operation is not differentiable. The loss vector needs to be resorted at every iteration because the tail can change with every shares update.

I even rewrote the TVar calculation function to avoid the sorting but the optimizer then fails to converge with reasonable constraints level, or it just ignores them.

If it helps, this is what I would normally do in Python using Scipy

constraints = [
{'type': 'ineq', 'fun': lambda x: tvar_down(x)},
{'type': 'ineq', 'fun': lambda x: tvar_up(x)},
]
2

Unfortunately you can’t formulate this constraint in a purely continuous linear program. You will have to either add binary variables with a formulation trick like optimization - Linear constraint expressing the sum of the $k$ largest elements - Mathematics Stack Exchange, or pick another modeling paradigm.

1 Like
3

Hi @Alessio_Bellisomi,

The paper you’re looking for is:

Rockafellar, R. Tyrrell, and Stanislav Uryasev. “Optimization of conditional value-at-risk.” Journal of risk 2 (2000): 21-42.

If you change rp from a number to the quantile 0.2, this is equivalent to CVaR.

julia> using JuMP, Ipopt, Statistics

julia> function min_cvar_fn(x, alpha)
           model = Model(Ipopt.Optimizer)
           set_silent(model)
           @variable(model, z[1:length(x)] >= 0)
           @variable(model, var)
           @constraint(model, z .>= x .- var)
           @objective(model, Min, var + 1 / alpha * sum(z) / length(z))
           optimize!(model)
           @assert is_solved_and_feasible(model)
           return objective_value(model)
       end
min_cvar_fn (generic function with 1 method)

julia> x = rand(10)
10-element Vector{Float64}:
 0.2735943053314943
 0.9101308567374451
 0.0874812908694732
 0.8289752673157755
 0.42123006421472475
 0.3416336067914879
 0.23804878901581594
 0.9742698820738032
 0.9221506558571922
 0.47265903686773114

julia> mean(sort(x; rev = true)[1:2])
0.9482102689654976

julia> min_cvar_fn(x, 0.2)
0.9482102415031161

You probably want something like:

model = Model(Ipopt.Optimizer)
set_silent(model)
@variable(model, minimum_shares[i] <= x[i in 1:num_contracts] <= maximum_shares[i])
set_start_value.(x,  shares)
@variable(model, portfolio_losses[1:num_scenarios])
@constraint(model, portfolio_losses .== losses_matrix' * x)
@variable(model, var)
@variable(model, tvar_lower_bound <= cvar <= tvar_upper_bound)
@variable(model, z[i in 1:num_scenarios] >= 0)
@constraint(model, z .>= portfolio_losses .- var)
@constraint(model, cvar == var + 1 / alpha * sum(z) / num_scenarios)
@objective(model, Max, objective_prem_vs_initial(x))
optimize!(model)
is_success = is_solved_and_feasible(model)
2 Likes
4

that is really interesting, thank you for helping out!

1 Like