function non_convex_opt(net_data)
    """
    solves non-convex gas network optimization
    """
    # build model
    model = Model(optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0))
    # variable declaration
    @variable(model, ğ�›‘[1:length(net_data[:N])])
    @variable(model, Ï†[1:length(net_data[:E])])
    @variable(model, Îº[1:length(net_data[:E])])
    @variable(model, Ï‘[1:length(net_data[:N])])
    # minimize gas production cost
    @objective(model, Min, Ï‘'*diagm(vec(net_data[:c]))*Ï‘)
    # gas variable limits
    @constraint(model, inj_lim_max[i=net_data[:N]], Ï‘[i] <= net_data[:Ï‘Ì…][i])
    @constraint(model, inj_lim_min[i=net_data[:N]], Ï‘[i] >= net_data[:Ï‘Ì²][i])
    @constraint(model, pre_lim_max[i=net_data[:N]], ğ�›‘[i] <= net_data[:Ï�Ì…][i]^2)
    @constraint(model, pre_lim_min[i=net_data[:N]], ğ�›‘[i] >= net_data[:Ï�Ì²][i]^2)
    @constraint(model, com_lim_max[i=net_data[:E]], Îº[i] <= net_data[:ÎºÌ…][i])
    @constraint(model, com_lim_min[i=net_data[:E]], Îº[i] >= net_data[:ÎºÌ²][i])
    # gas flow equations
    @NLconstraint(model, w_eq[l=net_data[:E]], Ï†[l]*abs(Ï†[l]) - net_data[:k][l]^2 *(ğ�›‘[ns(l)] + Îº[l] - ğ�›‘[nr(l)]) == 0)
    @constraint(model, gas_bal, Ï‘ .- net_data[:Î´] .- net_data[:B]*Îº .== net_data[:A]*Ï†)
    @constraint(model, Ï†_pl[l=net_data[:E_a]], Ï†[l] >= 0)
    # solve model
    optimize!(model)
    @info("non_convex model terminates with status: $(termination_status(model))")
    # return solution
    solution = Dict(
    :Ï‘ => JuMP.value.(Ï‘),
    :ğ�›‘ => JuMP.value.(ğ�›‘),
    :Ï† => JuMP.value.(Ï†),
    :Îº => JuMP.value.(Îº),
    :cost => JuMP.objective_value.(model),
    :model => model)
    return solution
end

function linearized_opt(net_data,lin_res)
    """
    solves linearized gas network optimization
    """
    # build model
    model = Model(optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0))
    # variable declaration
    @variable(model, ğ�›‘[1:length(net_data[:N])])
    @variable(model, Ï†[1:length(net_data[:E])])
    @variable(model, Îº[1:length(net_data[:E])])
    @variable(model, Ï‘[1:length(net_data[:N])])
    # minimize gas production cost
    @objective(model, Min, Ï‘'*diagm(vec(net_data[:c]))*Ï‘)
    # gas variable limits
    @constraint(model, inj_lim_max[i=net_data[:N]], Ï‘[i] <= net_data[:Ï‘Ì…][i])
    @constraint(model, inj_lim_min[i=net_data[:N]], Ï‘[i] >= net_data[:Ï‘Ì²][i])
    @constraint(model, pre_lim_max[i=net_data[:N]], ğ�›‘[i] <= net_data[:Ï�Ì…][i]^2)
    @constraint(model, pre_lim_min[i=net_data[:N]], ğ�›‘[i] >= net_data[:Ï�Ì²][i]^2)
    @constraint(model, com_lim_max[i=net_data[:E]], Îº[i] <= net_data[:ÎºÌ…][i])
    @constraint(model, com_lim_min[i=net_data[:E]], Îº[i] >= net_data[:ÎºÌ²][i])
    # gas flow equations
    @constraint(model, w_eq, Ï† .== lin_res[:Ï‚1] + lin_res[:Ï‚2] * ğ�›‘ + lin_res[:Ï‚3] * Îº)
    @constraint(model, gas_bal, Ï‘ .- net_data[:Î´] .- net_data[:B]*Îº .== net_data[:A]*Ï†)
    @constraint(model, Ï†_pl[l=net_data[:E_a]], Ï†[l] >= 0)
    @constraint(model, ğ�›‘[net_data[:ref]] == lin_res[:ğ�›‘Ì‡][net_data[:ref]])
    # solve model
    optimize!(model)
    @info("linearized model terminates with status: $(termination_status(model))")
    # return solution
    solution = Dict(
    :Ï‘ => JuMP.value.(Ï‘),
    :ğ�›‘ => JuMP.value.(ğ�›‘),
    :Ï† => JuMP.value.(Ï†),
    :Îº => JuMP.value.(Îº),
    :cost => JuMP.objective_value.(model),
    :Î»_c => JuMP.dual.(gas_bal),
    :Î»_w => JuMP.dual.(w_eq),
    :model => model)
    return solution
end

function chance_con_opt(net_data,lin_res,forecast,settings)
    """
    solves chance-constrained gas network optimization
    """
    # build model
    model = Model(optimizer_with_attributes(Mosek.Optimizer, "LOG" => 0))
    # variable declaration
    @variable(model, ğ�›‘[1:length(net_data[:N])])
    @variable(model, Ï†[1:length(net_data[:E])])
    @variable(model, Îº[1:length(net_data[:E])])
    @variable(model, Ï‘[1:length(net_data[:N])])
    @variable(model, Î±[1:length(net_data[:N]),1:length(net_data[:N])])
    @variable(model, Î²[1:length(net_data[:E]),1:length(net_data[:N])])
    @variable(model, s_ğ�›‘[1:length(net_data[:N])])
    @variable(model, s_Ï†[1:length(net_data[:E])])
    @variable(model, c_Ï‘[1:length(net_data[:N])])
    @variable(model, c_Î±[1:length(net_data[:N])])
    # minimize expected gas injection cost + variance penalty
    @objective(model, Min, sum(c_Ï‘) + sum(c_Î±) + settings[:Ïˆ_ğ�›‘]*sum(s_ğ�›‘) + settings[:Ïˆ_Ï†]*sum(s_Ï†))
    # exp cost quadratic terms
    @constraint(model, Î»_Î¼_u_Ï‘[n=net_data[:N]], [1/2;c_Ï‘[n];net_data[:cÌ€][n]*Ï‘[n]] in RotatedSecondOrderCone())
    @constraint(model, Î»_Î¼_u_Î±[n=net_data[:N]], [1/2;c_Î±[n];net_data[:cÌ€][n]*forecast[:Î£Â½]*Î±[n,:]] in RotatedSecondOrderCone())
    # gas flow equations
    @constraint(model, Î»_c, 0 .== net_data[:A]*Ï† .- Ï‘ .+ net_data[:B]*Îº .+ net_data[:Î´])
    @constraint(model, Î»_w, 0 .== Ï† .- lin_res[:Ï‚1] .- lin_res[:Ï‚2] * ğ�›‘ .- lin_res[:Ï‚3] * Îº)
    @constraint(model, Î»_Ï€, ğ�›‘[net_data[:ref]] == lin_res[:ğ�›‘Ì‡][net_data[:ref]])
    @constraint(model, Î»_r, (Î± - net_data[:B]*Î²)'*ones(length(net_data[:N])) .== forecast[:I_Î´])
    # # chance constraints
    # variance control
    @constraint(model, Î»_u_s_Ï€[n=net_data[:N]],  [s_ğ�›‘[n] - 0; forecast[:Î£Â½] * (lin_res[:Ï‚Ì†2]*(Î± .- lin_res[:Ï‚Ì‚3]*Î² .- diagm(ones(length(net_data[:N])))))[n,:]] in SecondOrderCone())
    @constraint(model, Î»_u_s_Ï†[l=net_data[:E]],  [s_Ï†[l] - 0; forecast[:Î£Â½] * (lin_res[:Ï‚Ì€2]*(Î± .- diagm(ones(length(net_data[:N])))) - lin_res[:Ï‚Ì€3]*Î²)[l,:]] in SecondOrderCone())
    # pressure limits
    @constraint(model, Î»_u_Ï€Ì…[n=net_data[:N]],    [net_data[:Ï�Ì…][n]^2 - ğ�›‘[n]; Î¦(settings,net_data) * forecast[:Î£Â½] * (lin_res[:Ï‚Ì†2]*(Î± .- lin_res[:Ï‚Ì‚3]*Î² .- diagm(ones(length(net_data[:N])))))[n,:]] in SecondOrderCone())
    @constraint(model, Î»_u_Ï€Ì²[n=net_data[:N]],    [ğ�›‘[n] - net_data[:Ï�Ì²][n]^2; Î¦(settings,net_data) * forecast[:Î£Â½] * (lin_res[:Ï‚Ì†2]*(Î± .- lin_res[:Ï‚Ì‚3]*Î² .- diagm(ones(length(net_data[:N])))))[n,:]] in SecondOrderCone())
    # flow limits
    @constraint(model, Î»_u_Ï†Ì²[l=net_data[:E_a]],  [Ï†[l] - 0; Î¦(settings,net_data) * forecast[:Î£Â½] * (lin_res[:Ï‚Ì€2]*(Î± .- diagm(ones(length(net_data[:N])))) - lin_res[:Ï‚Ì€3]*Î²)[l,:]] in SecondOrderCone())
    # injection limits
    @constraint(model, Î»_u_Ï‘Ì…[n=net_data[:N]],  [net_data[:Ï‘Ì…][n] - Ï‘[n]; Î¦(settings,net_data) * forecast[:Î£Â½] * Î±[n,:]] in SecondOrderCone())
    @constraint(model, Î»_u_Ï‘Ì²[n=net_data[:N]],  [Ï‘[n] - net_data[:Ï‘Ì²][n]; Î¦(settings,net_data) * forecast[:Î£Â½] * Î±[n,:]] in SecondOrderCone())
    # compression limits
    @constraint(model, Î»_u_ÎºÌ…[l=net_data[:E]],  [net_data[:ÎºÌ…][l] - Îº[l]; Î¦(settings,net_data) * forecast[:Î£Â½] * Î²[l,:]] in SecondOrderCone())
    @constraint(model, Î»_u_ÎºÌ²[l=net_data[:E]],  [Îº[l] - net_data[:ÎºÌ²][l]; Î¦(settings,net_data) * forecast[:Î£Â½] * Î²[l,:]] in SecondOrderCone())
    # aux constraints
    @constraint(model, Î´Ì¸_Î², Î²[:,forecast[:N_Î´Ì¸]] .== 0)
    @constraint(model, Î´Ì¸_Î±, Î±[:,forecast[:N_Î´Ì¸]] .== 0)
    settings[:comp] == false ? @constraint(model, Î²[findall(x->x>0, net_data[:ÎºÌ…]),:] .==0) : NaN    # zero compressor response to uncertainty
    settings[:valv] == false ? @constraint(model, Î²[findall(x->x<0, net_data[:ÎºÌ²]),:] .==0) : NaN    # zero valve response to uncertainty
    # solve model
    optimize!(model)
    @info("policy optimization terminates with status: $(termination_status(model))")
    # return solution
    solution = Dict(
    # primal solution
    :obj => JuMP.objective_value.(model),
    :c_Ï‘ => JuMP.value.(c_Ï‘),
    :c_Î± => JuMP.value.(c_Î±),
    :Ï‘ => JuMP.value.(Ï‘),
    :ğ�›‘ => JuMP.value.(ğ�›‘),
    :Ï† => JuMP.value.(Ï†),
    :Îº => JuMP.value.(Îº),
    :Î± => JuMP.value.(Î±),
    :Î² => JuMP.value.(Î²),
    :s_ğ�›‘ => JuMP.value.(s_ğ�›‘),
    :s_Ï† => JuMP.value.(s_Ï†),
    :cost => JuMP.value.(Ï‘)'*diagm(vec(net_data[:c]))*JuMP.value.(Ï‘) + tr(JuMP.value.(Î±)'diagm(vec(net_data[:c]))JuMP.value.(Î±)*forecast[:Î£]),
    :Î»_Î¼_u_Ï‘ => JuMP.dual.(Î»_Î¼_u_Ï‘),
    :Î»_Î¼_u_Î± => JuMP.dual.(Î»_Î¼_u_Î±),
    :Î»_c => JuMP.dual.(Î»_c),
    :Î»_w => JuMP.dual.(Î»_w),
    :Î»_Ï€ => JuMP.dual.(Î»_Ï€),
    :Î»_r => JuMP.dual.(Î»_r),
    :Î»_u_Ï€Ì… => JuMP.dual.(Î»_u_Ï€Ì…),
    :Î»_u_Ï€Ì² => JuMP.dual.(Î»_u_Ï€Ì²),
    :Î»_u_Ï†Ì² => JuMP.dual.(Î»_u_Ï†Ì²),
    :Î»_u_s_Ï€ => JuMP.dual.(Î»_u_s_Ï€),
    :Î»_u_s_Ï† => JuMP.dual.(Î»_u_s_Ï†),
    :Î»_u_Ï‘Ì… => JuMP.dual.(Î»_u_Ï‘Ì…),
    :Î»_u_Ï‘Ì² => JuMP.dual.(Î»_u_Ï‘Ì²),
    :Î»_u_ÎºÌ… => JuMP.dual.(Î»_u_ÎºÌ…),
    :Î»_u_ÎºÌ² => JuMP.dual.(Î»_u_ÎºÌ²),
    :model => model)
    return solution
end

function stochastic_dual_solution(net_data,sol_stochastic,lin_res,forecast,settings)
    """
    extracts and processes the dual solution of the chance-constrained gas network optimization
    """
    # extract dual variables
    Î»_c = sol_stochastic[:Î»_c]
    Î»_w = sol_stochastic[:Î»_w]
    Î»_Ï€Ì‡ = sol_stochastic[:Î»_Ï€]
    Î»_r = sol_stochastic[:Î»_r]
    Î»_Ï‘ = [sol_stochastic[:Î»_Î¼_u_Ï‘][i][1] for i in net_data[:N]]
    Î»_Î± = [sol_stochastic[:Î»_Î¼_u_Î±][i][1] for i in net_data[:N]]
    Î»_Ï€Ì… = [sol_stochastic[:Î»_u_Ï€Ì…][i][1] for i in net_data[:N]]
    Î»_Ï€Ì² = [sol_stochastic[:Î»_u_Ï€Ì²][i][1] for i in net_data[:N]]
    Î»_Ï‘Ì… = [sol_stochastic[:Î»_u_Ï‘Ì…][i][1] for i in net_data[:N]]
    Î»_Ï‘Ì² = [sol_stochastic[:Î»_u_Ï‘Ì²][i][1] for i in net_data[:N]]
    Î»_ÎºÌ… = [sol_stochastic[:Î»_u_ÎºÌ…][i][1] for i in net_data[:E]]
    Î»_ÎºÌ² = [sol_stochastic[:Î»_u_ÎºÌ²][i][1] for i in net_data[:E]]
    Î»_s_Ï€ = [sol_stochastic[:Î»_u_s_Ï€][i][1] for i in net_data[:N]]
    Î»_s_Ï† = [sol_stochastic[:Î»_u_s_Ï†][i][1] for i in net_data[:E]]
    u_Ï‘ = [sol_stochastic[:Î»_Î¼_u_Ï‘][i][3] for i in net_data[:N]]
    u_Î± = [sol_stochastic[:Î»_Î¼_u_Î±][i][j] for i in net_data[:N], j in 3:length(net_data[:N])+2]
    u_Ï€Ì… = [sol_stochastic[:Î»_u_Ï€Ì…][i][j] for i in net_data[:N], j in 2:length(net_data[:N])+1]
    u_Ï€Ì² = [sol_stochastic[:Î»_u_Ï€Ì²][i][j] for i in net_data[:N], j in 2:length(net_data[:N])+1]
    u_Ï‘Ì… = [sol_stochastic[:Î»_u_Ï‘Ì…][i][j] for i in net_data[:N], j in 2:length(net_data[:N])+1]
    u_Ï‘Ì² = [sol_stochastic[:Î»_u_Ï‘Ì²][i][j] for i in net_data[:N], j in 2:length(net_data[:N])+1]
    u_ÎºÌ… = [sol_stochastic[:Î»_u_ÎºÌ…][i][j] for i in net_data[:E], j in 2:length(net_data[:N])+1]
    u_ÎºÌ² = [sol_stochastic[:Î»_u_ÎºÌ²][i][j] for i in net_data[:E], j in 2:length(net_data[:N])+1]
    u_s_Ï€ = [sol_stochastic[:Î»_u_s_Ï€][i][j] for i in net_data[:N], j in 2:length(net_data[:N])+1]
    u_s_Ï† = [sol_stochastic[:Î»_u_s_Ï†][i][j] for i in net_data[:E], j in 2:length(net_data[:N])+1]
    Î»_Ï†Ì² = zeros(length(net_data[:E]))
    u_Ï†Ì² = zeros(length(net_data[:E]),length(net_data[:N]))
    for i in net_data[:E]
        i âˆˆ net_data[:E_a] ? Î»_Ï†Ì²[i] = sol_stochastic[:Î»_u_Ï†Ì²][i][1] : NaN
        i âˆˆ net_data[:E_a] ? u_Ï†Ì²[i,:] = sol_stochastic[:Î»_u_Ï†Ì²][i][2:end] : NaN
    end
    Î¼_Ï‘ = [sol_stochastic[:Î»_Î¼_u_Ï‘][i][2] for i in net_data[:N]]
    Î¼_Î± =  [sol_stochastic[:Î»_Î¼_u_Î±][i][2] for i in net_data[:N]]

    # partial lagrangian function
    L_part = (
    + Î»_c'*(net_data[:A]*sol_stochastic[:Ï†] .- sol_stochastic[:Ï‘] + net_data[:B]*sol_stochastic[:Îº] + net_data[:Î´])
    + Î»_r'*(ones(length(net_data[:N])) - (sol_stochastic[:Î±]-net_data[:B]*sol_stochastic[:Î²])'*ones(length(net_data[:N])))
    + Î»_w'*(sol_stochastic[:Ï†] - lin_res[:Ï‚1] - lin_res[:Ï‚2]*sol_stochastic[:ğ�›‘] - lin_res[:Ï‚3]*sol_stochastic[:Îº])
    - Î»_s_Ï†'*sol_stochastic[:s_Ï†]
    - Î»_s_Ï€'*sol_stochastic[:s_ğ�›‘]
    - Î»_Ï†Ì²'*sol_stochastic[:Ï†]
    - Î»_Ï€Ì…'*(net_data[:Ï�Ì…].^2 - sol_stochastic[:ğ�›‘])
    - Î»_Ï€Ì²'*(sol_stochastic[:ğ�›‘] - net_data[:Ï�Ì²].^2)
    - sum((u_s_Ï† + Î¦(settings,net_data)*u_Ï†Ì²)[l,:]'*forecast[:Î£Â½]*(lin_res[:Ï‚Ì€2]*(sol_stochastic[:Î±] - diagm(ones(length(net_data[:N])))) - lin_res[:Ï‚Ì€3]*sol_stochastic[:Î²])[l,:] for l in net_data[:E])
    - sum((u_s_Ï€ + Î¦(settings,net_data)*u_Ï€Ì… + Î¦(settings,net_data)*u_Ï€Ì²)[n,:]' * forecast[:Î£Â½]*(lin_res[:Ï‚Ì†2]*(sol_stochastic[:Î±] - lin_res[:Ï‚Ì‚3]*sol_stochastic[:Î²] - diagm(ones(length(net_data[:N])))))[n,:] for n in net_data[:N])
    )

    # partial lagrangian function decomposition
    R_free = + Î»_w' * lin_res[:Ï‚1]
    R_rent = (
            - Î»_c'*net_data[:A] * sol_stochastic[:Ï†]
            - Î»_w' * sol_stochastic[:Ï†]
            + Î»_w' * lin_res[:Ï‚2] * sol_stochastic[:ğ�›‘]
            + Î»_Ï†Ì²'*sol_stochastic[:Ï†]
            + Î»_Ï€Ì…'*(net_data[:Ï�Ì…].^2 - sol_stochastic[:ğ�›‘])
            + Î»_Ï€Ì²'*(sol_stochastic[:ğ�›‘] - net_data[:Ï�Ì²].^2)
            + Î»_s_Ï†'*sol_stochastic[:s_Ï†]
            + Î»_s_Ï€'*sol_stochastic[:s_ğ�›‘]
            )
    R_inj = (
            + Î»_c' * sol_stochastic[:Ï‘]
            + Î»_r' * sol_stochastic[:Î±]'*ones(length(net_data[:N]))
            + sum((u_s_Ï† + Î¦(settings,net_data)*u_Ï†Ì²)[l,:]'*forecast[:Î£Â½]*(lin_res[:Ï‚Ì€2]*sol_stochastic[:Î±])[l,:] for l in net_data[:E])
            + sum((u_s_Ï€ + Î¦(settings,net_data)*u_Ï€Ì… + Î¦(settings,net_data)*u_Ï€Ì²)[n,:]' * forecast[:Î£Â½]*(lin_res[:Ï‚Ì†2]*sol_stochastic[:Î±])[n,:] for n in net_data[:N])
            )
    R_act = (
            - Î»_c' * net_data[:B] * sol_stochastic[:Îº]
            - Î»_r' * (net_data[:B]*sol_stochastic[:Î²])' * ones(length(net_data[:N]))
            + Î»_w' * lin_res[:Ï‚3] * sol_stochastic[:Îº]
            - sum((u_s_Ï† + Î¦(settings,net_data)*u_Ï†Ì²)[l,:]'*forecast[:Î£Â½]*(lin_res[:Ï‚Ì€3]*sol_stochastic[:Î²])[l,:] for l in net_data[:E])
            - sum((u_s_Ï€ + Î¦(settings,net_data)*u_Ï€Ì… + Î¦(settings,net_data)*u_Ï€Ì²)[n,:]' * forecast[:Î£Â½]*(lin_res[:Ï‚Ì†2]*lin_res[:Ï‚Ì‚3]*sol_stochastic[:Î²])[n,:] for n in net_data[:N])
            )
    R_con = (
            + Î»_c' * net_data[:Î´]
            + Î»_r' * ones(length(net_data[:N]))
            + sum((u_s_Ï† + Î¦(settings,net_data)*u_Ï†Ì²)[l,:]'*forecast[:Î£Â½]*lin_res[:Ï‚Ì€2][l,:] for l in net_data[:E])
            + sum((u_s_Ï€ + Î¦(settings,net_data)*u_Ï€Ì… + Î¦(settings,net_data)*u_Ï€Ì²)[n,:]' * forecast[:Î£Â½]*lin_res[:Ï‚Ì†2][n,:] for n in net_data[:N])
            )

    # revenue balance
    R_con - R_inj - R_act - R_rent - R_free
    # revenue decomposition
    R_inj_nom_bal = Î»_c' * sol_stochastic[:Ï‘]
    R_inj_rec_bal = Î»_r' * sol_stochastic[:Î±]'*ones(length(net_data[:N]))
    R_inj_net_lim = sum(Î¦(settings,net_data)*u_Ï†Ì²[l,:]'*forecast[:Î£Â½]*(lin_res[:Ï‚Ì€2]*sol_stochastic[:Î±])[l,:] for l in net_data[:E]) + sum((Î¦(settings,net_data)*u_Ï€Ì… + Î¦(settings,net_data)*u_Ï€Ì²)[n,:]' * forecast[:Î£Â½]*(lin_res[:Ï‚Ì†2]*sol_stochastic[:Î±])[n,:] for n in net_data[:N])
    R_inj_net_var = sum(u_s_Ï†[l,:]'*forecast[:Î£Â½]*(lin_res[:Ï‚Ì€2]*sol_stochastic[:Î±])[l,:] for l in net_data[:E]) + sum(u_s_Ï€[n,:]' * forecast[:Î£Â½]*(lin_res[:Ï‚Ì†2]*sol_stochastic[:Î±])[n,:] for n in net_data[:N])

    R_act_nom_bal = Î»_w' * lin_res[:Ï‚3] * sol_stochastic[:Îº] - Î»_c' * net_data[:B] * sol_stochastic[:Îº]
    R_act_rec_bal = - Î»_r' * (net_data[:B]*sol_stochastic[:Î²])' * ones(length(net_data[:N]))
    R_act_net_lim = - sum((Î¦(settings,net_data)*u_Ï†Ì²)[l,:]'*forecast[:Î£Â½]*(lin_res[:Ï‚Ì€3]*sol_stochastic[:Î²])[l,:] for l in net_data[:E]) - sum((Î¦(settings,net_data)*u_Ï€Ì… + Î¦(settings,net_data)*u_Ï€Ì²)[n,:]' * forecast[:Î£Â½]*(lin_res[:Ï‚Ì†2]*lin_res[:Ï‚Ì‚3]*sol_stochastic[:Î²])[n,:] for n in net_data[:N])
    R_act_net_var = - sum(u_s_Ï†[l,:]'*forecast[:Î£Â½]*(lin_res[:Ï‚Ì€3]*sol_stochastic[:Î²])[l,:] for l in net_data[:E]) - sum(u_s_Ï€[n,:]' * forecast[:Î£Â½]*(lin_res[:Ï‚Ì†2]*lin_res[:Ï‚Ì‚3]*sol_stochastic[:Î²])[n,:] for n in net_data[:N])

    R_con_nom_bal = Î»_c' * net_data[:Î´]
    R_con_rec_bal = Î»_r' * ones(length(net_data[:N]))
    R_con_net_lim = sum((Î¦(settings,net_data)*u_Ï†Ì²)[l,:]'*forecast[:Î£Â½]*lin_res[:Ï‚Ì€2][l,:] for l in net_data[:E]) + sum((Î¦(settings,net_data)*u_Ï€Ì… + Î¦(settings,net_data)*u_Ï€Ì²)[n,:]' * forecast[:Î£Â½]*lin_res[:Ï‚Ì†2][n,:] for n in net_data[:N])
    R_con_net_var = sum(u_s_Ï†[l,:]'*forecast[:Î£Â½]*lin_res[:Ï‚Ì€2][l,:] for l in net_data[:E]) + sum(u_s_Ï€[n,:]' * forecast[:Î£Â½]*lin_res[:Ï‚Ì†2][n,:] for n in net_data[:N])

    Revenue_decomposition = DataFrame(source = ["nom_bal","rec_bal","net_lim","net_var","total"],
                                            inj = [R_inj_nom_bal,R_inj_rec_bal,R_inj_net_lim,R_inj_net_var, R_inj_nom_bal+R_inj_rec_bal+R_inj_net_lim+R_inj_net_var],
                                            act = [R_act_nom_bal,R_act_rec_bal,R_act_net_lim,R_act_net_var, R_act_nom_bal+R_act_rec_bal+R_act_net_lim+R_act_net_var],
                                            con = [R_con_nom_bal,R_con_rec_bal,R_con_net_lim,R_con_net_var, R_con_nom_bal+R_con_rec_bal+R_con_net_lim+R_con_net_var]
    )

    # compute indvidual revenues & profits
    R_ind_inj = [Î»_c[n]*sol_stochastic[:Ï‘][n] + (Î»_r' + lin_res[:Ï‚Ì†2][:,n]'*(u_s_Ï€ .+ Î¦(settings,net_data) * u_Ï€Ì… .+ Î¦(settings,net_data)*  u_Ï€Ì²)*forecast[:Î£Â½] + lin_res[:Ï‚Ì€2][:,n]'*(u_s_Ï† .+ Î¦(settings,net_data) * u_Ï†Ì²)*forecast[:Î£Â½])*sol_stochastic[:Î±][n,:] for n in net_data[:N]]
    Î _inj = [R_ind_inj[n] - sol_stochastic[:c_Ï‘][n] - sol_stochastic[:c_Î±][n] for n in net_data[:N]]
    R_ind_act = [ -Î»_c'*net_data[:B][:,l] * sol_stochastic[:Îº][l] + lin_res[:Ï‚3][:,l]' * Î»_w * sol_stochastic[:Îº][l] - ones(length(net_data[:N]))' * net_data[:B][:,l] * Î»_r' * sol_stochastic[:Î²][l,:] - (lin_res[:Ï‚Ì†2]*lin_res[:Ï‚Ì‚3])[:,l]' * (u_s_Ï€ .+ Î¦(settings,net_data) * u_Ï€Ì… .+ Î¦(settings,net_data)*  u_Ï€Ì²) * forecast[:Î£Â½] * sol_stochastic[:Î²][l,:] - lin_res[:Ï‚Ì€3][:,l]' * (u_s_Ï† .+ Î¦(settings,net_data) * u_Ï†Ì²) * forecast[:Î£Â½] * sol_stochastic[:Î²][l,:] for l in net_data[:E_a]]
    R_ind_rent = (- Î»_c'*net_data[:A] * sol_stochastic[:Ï†] - Î»_w' * sol_stochastic[:Ï†] + Î»_w' * lin_res[:Ï‚2] * sol_stochastic[:ğ�›‘] + Î»_Ï†Ì²'*sol_stochastic[:Ï†] + Î»_Ï€Ì…'*(net_data[:Ï�Ì…].^2 - sol_stochastic[:ğ�›‘]) + Î»_Ï€Ì²'*(sol_stochastic[:ğ�›‘] - net_data[:Ï�Ì²].^2) + Î»_s_Ï†'*sol_stochastic[:s_Ï†] + Î»_s_Ï€'*sol_stochastic[:s_ğ�›‘])
    # R_ind_con = [Î»_c[n]*net_data[:Î´][n] + Î»_r[n] + sum(lin_res[:Ï‚Ì€2][l,n]*(u_s_Ï† .+ Î¦(settings,net_data) * u_Ï†Ì²)[l,:]'*forecast[:Î£Â½][n,:] for l in net_data[:E]) + sum(lin_res[:Ï‚Ì†2][k,n]*(u_s_Ï€ .+ Î¦(settings,net_data) * u_Ï€Ì… .+ Î¦(settings,net_data)*  u_Ï€Ì²)[k,:]'*forecast[:Î£Â½][n,:] for k in net_data[:N]) for n in net_data[:N]]
    R_ind_con = [Î»_c[n]*net_data[:Î´][n] + Î»_r[n] + forecast[:Î£Â½][n,:]'*(u_s_Ï† .+ Î¦(settings,net_data) * u_Ï†Ì²)'*lin_res[:Ï‚Ì€2][:,n] + forecast[:Î£Â½][n,:]'*(u_s_Ï€ .+ Î¦(settings,net_data) * u_Ï€Ì… .+ Î¦(settings,net_data)*  u_Ï€Ì²)'*lin_res[:Ï‚Ì†2][:,n] for n in net_data[:N]]

    # check stationarity conditions
    âˆ‚L_âˆ‚s_Ï† = settings[:Ïˆ_Ï†] .- Î»_s_Ï†
    âˆ‚L_âˆ‚s_Ï€ = settings[:Ïˆ_ğ�›‘] .- Î»_s_Ï€
    âˆ‚L_âˆ‚c_Ï‘ = 1 .- Î¼_Ï‘
    âˆ‚L_âˆ‚c_Î± = 1 .- Î¼_Î±
    âˆ‚L_âˆ‚Ï€ = Î»_Ï€Ì… .- Î»_Ï€Ì² .- vec(Î»_w'*lin_res[:Ï‚2]) ; âˆ‚L_âˆ‚Ï€[net_data[:ref]] = âˆ‚L_âˆ‚Ï€[net_data[:ref]] - Î»_Ï€Ì‡
    âˆ‚L_âˆ‚Ï† = vec(Î»_c'*net_data[:A]) .+ Î»_w .- Î»_Ï†Ì²
    âˆ‚L_âˆ‚Ï‘ = - u_Ï‘ .* net_data[:cÌ€] .- Î»_c .+ Î»_Ï‘Ì… .- Î»_Ï‘Ì²
    âˆ‚L_âˆ‚Îº = vec(Î»_c'*net_data[:B]) .- vec(Î»_w'*lin_res[:Ï‚3]) .+ Î»_ÎºÌ… .- Î»_ÎºÌ²
    âˆ‚L_âˆ‚Î± = zeros(length(net_data[:N]),length(net_data[:N]))
    âˆ‚L_âˆ‚Î² = zeros(length(net_data[:E]),length(net_data[:N]))
    for n in net_data[:N]
        # âˆ‚L_âˆ‚Î±[n,:] = 2*net_data[:c][n]*forecast[:Î£]*sol_stochastic[:Î±][n,:] .- Î»_r .- forecast[:Î£Â½] * (u_Ï†' * lin_res[:Ï‚Ì€2][:,n] + u_Ï€' * lin_res[:Ï‚Ì†2][:,n] + Î¦(settings,net_data) * u_Ï‘Ì…[n,:] + Î¦(settings,net_data) * u_Ï‘Ì²[n,:])
        âˆ‚L_âˆ‚Î±[n,:] = -forecast[:Î£Â½]*u_Î±[n,:]*net_data[:cÌ€][n] .- Î»_r .- forecast[:Î£Â½] * ((u_s_Ï† .+ Î¦(settings,net_data) * u_Ï†Ì²)' * lin_res[:Ï‚Ì€2][:,n] + (u_s_Ï€ .+ Î¦(settings,net_data) * u_Ï€Ì… .+ Î¦(settings,net_data)*  u_Ï€Ì²)' * lin_res[:Ï‚Ì†2][:,n] + Î¦(settings,net_data) * u_Ï‘Ì…[n,:] + Î¦(settings,net_data) * u_Ï‘Ì²[n,:])
    end
    for l in net_data[:E]
        âˆ‚L_âˆ‚Î²[l,:] = ones(length(net_data[:N]))'*net_data[:B][:,l]*Î»_r .+ forecast[:Î£Â½] * (u_s_Ï† .+ Î¦(settings,net_data) * u_Ï†Ì²)' * lin_res[:Ï‚Ì€3][:,l] .+ forecast[:Î£Â½] * (u_s_Ï€ .+ Î¦(settings,net_data) * u_Ï€Ì… .+ Î¦(settings,net_data)*  u_Ï€Ì²)' * (lin_res[:Ï‚Ì†2]*lin_res[:Ï‚Ì‚3])[:,l] .- forecast[:Î£Â½] * Î¦(settings,net_data) * (u_ÎºÌ… + u_ÎºÌ²)[l,:]
    end
    # total stationarity conditions mismatch
    mismatch = (sum(abs.(âˆ‚L_âˆ‚c_Ï‘)) + sum(abs.(âˆ‚L_âˆ‚c_Î±)) + sum(abs.(âˆ‚L_âˆ‚s_Ï†))
                + sum(abs.(âˆ‚L_âˆ‚s_Ï€)) + sum(abs.(âˆ‚L_âˆ‚Ï€)) + sum(abs.(âˆ‚L_âˆ‚Ï†))
                + sum(abs.(âˆ‚L_âˆ‚Ï‘)) + sum(abs.(âˆ‚L_âˆ‚Îº)) + sum(abs.(âˆ‚L_âˆ‚Î²)))
    mismatch <= 1e-2 ? @info("stationarity conditions hold; total mismatch : $(round(mismatch,digits=8))") : @warn("stationarity conditions do not hold;  mismatch : $(round(mismatch,digits=8))")

    # compute dual objective function
    dual_obj = (Î»_c' * net_data[:Î´]
    + sum(Î»_r)
    - Î»_w' * lin_res[:Ï‚1]
    - 0.5 * sum(Î»_Ï‘)
    - 0.5 * sum(Î»_Î±)
    - Î»_Ï€Ì…' * net_data[:Ï�Ì…].^2
    + Î»_Ï€Ì²' * net_data[:Ï�Ì²].^2
    - Î»_Ï‘Ì…' * net_data[:Ï‘Ì…]
    + Î»_Ï‘Ì²' * net_data[:Ï‘Ì²]
    - Î»_ÎºÌ…' * net_data[:ÎºÌ…]
    + Î»_ÎºÌ²' * net_data[:ÎºÌ²]
    + sum((u_s_Ï† .+ Î¦(settings,net_data) * u_Ï†Ì²)[l,:]' * forecast[:Î£Â½] * (lin_res[:Ï‚Ì€2] * diagm(ones(length(net_data[:N]))))[l,:] for l in net_data[:E])
    + sum((u_s_Ï€ .+ Î¦(settings,net_data) * u_Ï€Ì… .+ Î¦(settings,net_data)*  u_Ï€Ì²)[n,:]' * forecast[:Î£Â½] * (lin_res[:Ï‚Ì†2] * diagm(ones(length(net_data[:N]))))[n,:] for n in net_data[:N])
    + Î»_Ï€Ì‡*lin_res[:ğ�›‘Ì‡][net_data[:ref]])
    # check duality gap
    duality_gap = norm(dual_obj-sol_stochastic[:obj])/sol_stochastic[:obj]*100
    duality_gap <= 1e-3 ? @info("strong duality holds; duality gap : $(round(duality_gap,digits=3))%") : @info("strong duality does not hold; duality gap : $(round(duality_gap,digits=3))%")

    # return dual solution
    return Dict(:dual_obj => dual_obj, :R_inj => round.(R_inj), :Î _inj => round.(Î _inj), :R_act => round.(R_act), :R_con => round.(R_con), :R_rent => round.(R_rent), :Revenue_decomposition => Revenue_decomposition)
end
