Orbits/julia_code/attempt2/SimulationSteps.jl

using LinearAlgebra
using Distributions
using StatsFuns
using Plots

abstract type  DestructionProbability end

struct LinearProbability <: DestructionProbability
    spontaneous
    intra_constellation
    inter_constellation
    debris
end
function (lp::LinearProbability)(Sᵢ, S, D,) 
    min(1,max(0, 
        lp.spontaneous + lp.intra_constellation * Sᵢ + lp.inter_constellation * (S-Sᵢ) + lp.debris * D)
    )
end

function individual_loss(Sᵢ, S, D,n=1,return_prob=false)
    #= 
        This function simulates how many satellites will be lost 
        from a given constellation.

        TODO: Requires calibration of some sort
    =#
    β₀ = 3e-4 #probability of spontaneous destruction
    β₁ = 1e-6 #proability of collision with satellite in same constellation
    β₂ = 1e-5 #probability of collision with satellite in different constellation
    β₃ = 1e-4 #probability of collision with debris

    #linear probability (with constraints)
    p =  min(1,max(0,β₀ + β₁ * Sᵢ + β₂ * (S-Sᵢ) + β₃ * D))
    #p =  logistic(β₀ + β₁ * Sᵢ + β₂ * (S-Sᵢ) + β₃ * D)

    dist = Binomial(Sᵢ, p)
    draw = rand(dist)

    if return_prob
        return draw, p
    else
        return draw
    end
end


function debris_transition(D,S⃗, Y⃗, X⃗, δ=0.005, Γ=0.5)
    #δ: 1/δ year decay rate
    #Γ: proportion of the mass of the satellite that is released into orbit during launch

    #new debris released equals the debris lost to deorbit + debris generated by collisions + debris due to launches 
    new_debris = (1-δ) * D .+ sum(Y⃗) .+ Γ*sum(X⃗)
    return new_debris
end


function stocks_transition(S⃗, Y⃗, X⃗)
    return S⃗ .+ X⃗ .- Y⃗
end


function overall_losses(S⃗, D)
    #=
    Simulate the losses for each constellation
    =#
    
    individual_outcomes = [ individual_loss(s, sum(S⃗),D,1,true) for s in S⃗ ]

    losses = []
    prob_distruction = []
    for (loss,prob) in individual_outcomes
        append!(losses,loss)
        append!(prob_distruction,prob)
    end

    return losses,prob_distruction

end

# The demand function for the market

function inverse_demand_cournot(S⃗,A=100)
    A - sum(S⃗)
end

function profit_function(s,x,C_operation,C_launch,T_launch, T_operation,Price)
    #= 
        The profit function each participant recieves
        Includes place for taxes and costs.
    =#
    (Price-T_operation-C_operation)*s - (C_launch + T_launch)*x
end





#Get histogram of loss numbers per 10_000
histogram(
    [sum(individual_loss(33,67,500,10_000)) for i in 1:10000]
    ,bins=range(0,200,length=201)
    ,label="Satellites lost per 10,000"
    )



#=
Basic simulations
1. What happens when both participants always launch a single satellite?
2. What happens when a single participant always launches a single satellite?
3. What happens when both participants maintain the cournot equilibrium?
=#

#Specify policies

abstract type Policy end
#=
Represent any policy as a struct holding parameters that can be called with 
state variables to get the chosen policy.
=#
struct ConstantLaunchRatePolicy <: Policy
    #=
    Operator launches fixed number of satellites per turn
    =#
    rate_numerator
end
(p::ConstantLaunchRatePolicy)(s,S,D) = p.rate_numerator

struct SustainmentPolicy <: Policy
    #=
    Operator launches enough satellites to maintain a certain size
    =#
    level
end
(p::SustainmentPolicy)(s,S,D) = max(0,p.level - s)

struct simulation_result
    state_path
    profit_path
    destruction_events
    launch_choices
    destruction_probabilities
end


function simulate_system(policies,starting_state,n)
    #=
    Given a set of policies and a starting state, calculate n steps ahead.
    =#

    state = fill(starting_state,n)
    profit_path = []
    decay_events=[]
    launch_choices = []
    destruction_probabilities = []

    for i in 2:n
        current_state = state[i-1]
        S⃗ = current_state[2:end]
        D = current_state[1]

        #Calculate actions from policies
        X = [policy(s,sum(S⃗),D) for (s,policy) in zip(S⃗,policies)]
        push!(launch_choices, X)

        #Information generation
        price = inverse_demand_cournot(S⃗,10000)
        push!(profit_path,[profit_function(s,x,4,100,0, 0,price) for (s,x) in zip(S⃗,X) ])

        #calculate transition to next state
        #loss of satellites
        (Y⃗,prob_distruction) = overall_losses(S⃗, D)
        push!(decay_events,Y⃗)
        push!(destruction_probabilities, prob_distruction)

        #debris
        new_debris = debris_transition(D, S⃗, Y⃗, X,0.05, )

        #stocks
        new_stocks = stocks_transition(S⃗,Y⃗,X)

        #update vector of states
        next_state = [new_debris,new_stocks...]

        #update vector
        state[i] = next_state
        println(current_state,next_state)
    end

    k = length(starting_state)

    return simulation_result(
        [ [v[i] for v in state] for i in 1:k], #states
        [ [v[i] for v in profit_path] for i in 1:3 ], #profit paths
        [ [v[i] for v in decay_events ] for i in 1:3], #counts of satellite distructions
        [ [v[i] for v in launch_choices ] for i in 1:3], #counts of satellite launches
        [ [v[i] for v in destruction_probabilities] for i in 1:3]
    )
end

#setup state
n=2000
starting_state = [10.0,0,0,0]
policies = [SustainmentPolicy(x) for x in 25:5:35]

simulation = simulate_system(policies,starting_state,n)


plot(
    simulation.state_path[2:4]
    ,labels=["stocks 1" "stocks 2" "Stocks 3"]
)

plot(
    simulation.state_path[1]
    ,labels="debris"
)


plot(
    simulation.state_path[1][3:end] .- simulation.state_path[1][2:end-1]
    ,labels="debris growth"
)

plot(    simulation.state_path 
    ,labels=["debris" "stocks 1" "stocks 2" "Stocks 3"]
)

plot(simulation.profit_path)

bar(
    simulation.destruction_events
    ,layout = grid(3,1)
    ,legend=:none
)

plot(
    simulation.launch_choices
    ,layout = grid(3,1)
    ,legend=:none
)

plot(simulation.destruction_probabilities)
#= Notes
With the (maintain equilibrium) approach, the probability 
of distruction rises to 100% within a certain amount of time.
=#

bar(
    [simulation.launch_choices, -1 .* simulation.destruction_events]
    ,layout = grid(3,1)
    ,ylims = (-25,25)
    #,label = ["Launches 1" "Launches 2" "Launches 3" "Destruction 1" "Destruction 2" "Destruction 3"]
    , legend = :none
    , title = ["Operator 1" "Operator 2" "Operator 3"]
    #,yaxis = (-20:5:20)
)


#Many simulations
sims = [simulate_system(policies,starting_state,n) for x in 1:1000]



state_plots = []
for sim in sims
    push!(state_plots,plot(sim.state_path[2:end]))
end

state_plots[2]

for p in state_plots
    display(p)
end


plot(sim1.state_path[2:end], color=:skyblue, alpha=0.5)


arr = [sim.state_path[2:end] for sim in sims]
plot(arr, color=:blue, alpha=0.01, legend=:none)


debarr = [ sim.state_path[1] for sim in sims] 
plot(debarr, color=:blue, alpha=0.05, legend=:none)

profitarr = [sim.profit_path for sim in sims]
plot(profitarr, color=:blue, alpha=0.01, legend=:none)


# Plot out probability of distruction