export cordoba_strg_tnep
""
function cordoba_strg_tnep(data::Dict{String,Any}, model_type::Type, solver; ref_extensions = [_PMACDC.add_ref_dcgrid!, _PMACDC.add_candidate_dcgrid!, _FP.add_candidate_storage!, _PM.ref_add_on_off_va_bounds!, _PM.ref_add_ne_branch!], setting = s, kwargs...)
    s = setting
    return _PM.run_model(data, model_type, solver, cordoba_post_strg_tnep; ref_extensions = [_PMACDC.add_ref_dcgrid!, _PMACDC.add_candidate_dcgrid!, _FP.add_candidate_storage!, _PM.ref_add_on_off_va_bounds!, _PM.ref_add_ne_branch!], setting = s, kwargs...)
end

# for distribution models
""
function cordoba_strg_tnep(data::Dict{String,Any}, model_type::Type{T}, optimizer; kwargs...) where T <: _PM.AbstractBFModel
    return _PM.run_model(data, model_type, optimizer, cordoba_post_strg_tnep;
                         ref_extensions = [add_candidate_storage!, _PM.ref_add_on_off_va_bounds!, ref_add_ne_branch_allbranches!, ref_add_frb_branch!, ref_add_oltc_branch!],
                         solution_processors = [_PM.sol_data_model!],
                         kwargs...)
end


# Here the problem is defined, which is then sent to the solver.
# It is basically a declarion of variables and constraint of the problem
""
function cordoba_post_strg_tnep(pm::_PM.AbstractPowerModel)
# VARIABLES: defined within PowerModels(ACDC) can directly be used, other variables need to be defined in the according sections of the code: see storage.jl
    for n in _PM.nw_ids(pm)
        _PM.variable_bus_voltage(pm; nw = n)
        _PM.variable_gen_power(pm; nw = n)
        _PM.variable_branch_power(pm; nw = n)
        _PM.variable_storage_power(pm; nw = n)

        _PMACDC.variable_voltage_slack(pm; nw = n)
        _PMACDC.variable_active_dcbranch_flow(pm; nw = n)
        _PMACDC.variable_dc_converter(pm; nw = n)
        _PMACDC.variable_dcbranch_current(pm; nw = n)
        _PMACDC.variable_dcgrid_voltage_magnitude(pm; nw = n)
        _FP.variable_absorbed_energy(pm; nw = n)
        _FP.variable_absorbed_energy_ne(pm; nw = n)


        # new variables for TNEP problem
        _PM.variable_ne_branch_indicator(pm; nw = n)
        _PM.variable_ne_branch_power(pm; nw = n)
        _PM.variable_ne_branch_voltage(pm; nw = n)
        _FP.variable_storage_power_ne(pm; nw = n)
        _PMACDC.variable_active_dcbranch_flow_ne(pm; nw = n)
        _PMACDC.variable_branch_ne(pm; nw = n)
        _PMACDC.variable_dc_converter_ne(pm; nw = n)
        _PMACDC.variable_dcbranch_current_ne(pm; nw = n)
        _PMACDC.variable_dcgrid_voltage_magnitude_ne(pm; nw = n)
    end
#OBJECTIVE see objective.jl
    objective_min_cost_storage_cordoba(pm)
#CONSTRAINTS: defined within PowerModels(ACDC) can directly be used, other constraints need to be defined in the according sections of the code: storage.jl
    for n in _PM.nw_ids(pm)
        _PM.constraint_model_voltage(pm; nw = n)
        _PM.constraint_ne_model_voltage(pm; nw = n)
        _PMACDC.constraint_voltage_dc(pm; nw = n)
        _PMACDC.constraint_voltage_dc_ne(pm; nw = n)
        for i in _PM.ids(pm, n, :ref_buses)
            _PM.constraint_theta_ref(pm, i, nw = n)
        end

        for i in _PM.ids(pm, n, :bus)
            _FP.constraint_power_balance_acne_dcne_strg(pm, i; nw = n)
        end
        if haskey(pm.setting, "allow_line_replacement") && pm.setting["allow_line_replacement"] == true
            for i in _PM.ids(pm, n, :branch)
                _FP.constraint_ohms_yt_from_repl(pm, i; nw = n)
                _FP.constraint_ohms_yt_to_repl(pm, i; nw = n)
                _FP.constraint_voltage_angle_difference_repl(pm, i; nw = n)
                _FP.constraint_thermal_limit_from_repl(pm, i; nw = n)
                _FP.constraint_thermal_limit_to_repl(pm, i; nw = n)
            end
        else
            for i in _PM.ids(pm, n, :branch)
                _PM.constraint_ohms_yt_from(pm, i; nw = n)
                _PM.constraint_ohms_yt_to(pm, i; nw = n)
                _PM.constraint_voltage_angle_difference(pm, i; nw = n)
                _PM.constraint_thermal_limit_from(pm, i; nw = n)
                _PM.constraint_thermal_limit_to(pm, i; nw = n)
            end
        end
        for i in _PM.ids(pm, n, :ne_branch)
            _PM.constraint_ne_ohms_yt_from(pm, i; nw = n)
            _PM.constraint_ne_ohms_yt_to(pm, i; nw = n)
            _PM.constraint_ne_voltage_angle_difference(pm, i; nw = n)
            _PM.constraint_ne_thermal_limit_from(pm, i; nw = n)
            _PM.constraint_ne_thermal_limit_to(pm, i; nw = n)
            if n > 1
                _PMACDC.constraint_candidate_acbranches_mp(pm, n, i)
            end
        end

        for i in _PM.ids(pm, n, :busdc)
            _PMACDC.constraint_power_balance_dc_dcne(pm, i; nw = n)
        end
        for i in _PM.ids(pm, n, :busdc_ne)
            _PMACDC.constraint_power_balance_dcne_dcne(pm, i; nw = n)
        end

        for i in _PM.ids(pm, n, :branchdc)
            _PMACDC.constraint_ohms_dc_branch(pm, i; nw = n)
        end
        for i in _PM.ids(pm, n, :branchdc_ne)
            _PMACDC.constraint_ohms_dc_branch_ne(pm, i; nw = n)
            _PMACDC.constraint_branch_limit_on_off(pm, i; nw = n)
            if n > 1
                _PMACDC.constraint_candidate_dcbranches_mp(pm, n, i)
            end
        end

        for i in _PM.ids(pm, n, :convdc)
            _PMACDC.constraint_converter_losses(pm, i; nw = n)
            _PMACDC.constraint_converter_current(pm, i; nw = n)
            _PMACDC.constraint_conv_transformer(pm, i; nw = n)
            _PMACDC.constraint_conv_reactor(pm, i; nw = n)
            _PMACDC.constraint_conv_filter(pm, i; nw = n)
            if pm.ref[:nw][n][:convdc][i]["islcc"] == 1
                _PMACDC.constraint_conv_firing_angle(pm, i; nw = n)
            end
        end
        for i in _PM.ids(pm, n, :convdc_ne)
            _PMACDC.constraint_converter_losses_ne(pm, i; nw = n)
            _PMACDC.constraint_converter_current_ne(pm, i; nw = n)
            _PMACDC.constraint_converter_limit_on_off(pm, i; nw = n)
            if n > 1
                _PMACDC.constraint_candidate_converters_mp(pm, n, i)
            end
            _PMACDC.constraint_conv_transformer_ne(pm, i; nw = n)
            _PMACDC.constraint_conv_reactor_ne(pm, i; nw = n)
            _PMACDC.constraint_conv_filter_ne(pm, i; nw = n)
            if pm.ref[:nw][n][:convdc_ne][i]["islcc"] == 1
                _PMACDC.constraint_conv_firing_angle_ne(pm, i; nw = n)
            end
        end

        for i in _PM.ids(pm, :storage, nw=n)
            _FP.constraint_storage_excl_slack(pm, i, nw = n)
            _PM.constraint_storage_thermal_limit(pm, i, nw = n)
            _PM.constraint_storage_losses(pm, i, nw = n)
        end
        for i in _PM.ids(pm, :ne_storage, nw=n)
            _FP.constraint_storage_excl_slack_ne(pm, i, nw = n)
            _FP.constraint_storage_thermal_limit_ne(pm, i, nw = n)
            _FP.constraint_storage_losses_ne(pm, i, nw = n)
            _FP.constraint_storage_bounds_ne(pm, i, nw = n)
        end
    end

    network_ids = sort(collect(_PM.nw_ids(pm)))
    n_1 = network_ids[1]
    n_last = network_ids[end]

    for i in _PM.ids(pm, :storage, nw = n_1)
        _FP.constraint_storage_state(pm, i, nw = n_1)
        _FP.constraint_maximum_absorption(pm, i, nw = n_1)
    end

    for i in _PM.ids(pm, :ne_storage, nw = n_1)
        _FP.constraint_storage_state_ne(pm, i, nw = n_1)
        _FP.constraint_maximum_absorption_ne(pm, i, nw = n_1)
    end

    for i in _PM.ids(pm, :storage, nw = n_last)
        _FP.constraint_storage_state_final(pm, i, nw = n_last)
    end

    for i in _PM.ids(pm, :ne_storage, nw = n_last)
        _FP.constraint_storage_state_final_ne(pm, i, nw = n_last)
    end

    for n_2 in network_ids[2:end]
        for i in _PM.ids(pm, :storage, nw = n_2)
            _FP.constraint_storage_state(pm, i, n_1, n_2)
            _FP.constraint_maximum_absorption(pm, i, n_1, n_2)
        end
        for i in _PM.ids(pm, :ne_storage, nw = n_2)
            _FP.constraint_storage_state_ne(pm, i, n_1, n_2)
            _FP.constraint_maximum_absorption_ne(pm, i, n_1, n_2)
            _FP.constraint_storage_investment(pm, n_1, n_2, i)
        end
        n_1 = n_2
    end

end
