diff --git a/julia_code/attempt2/SimulationSteps.jl b/julia_code/attempt2/SimulationSteps.jl index a3f8fcf..9efd496 100644 --- a/julia_code/attempt2/SimulationSteps.jl +++ b/julia_code/attempt2/SimulationSteps.jl @@ -3,27 +3,48 @@ using Distributions using StatsFuns using Plots -function individual_loss(Sᵢ, S, D,n=1) +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-5 #probability of spontaneous destruction + β₀ = 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 - p = β₀ + β₁ * Sᵢ + β₂ * (S-Sᵢ) + β₃ * D + #linear probability (with constraints) + p = min(1,max(0,β₀ + β₁ * Sᵢ + β₂ * (S-Sᵢ) + β₃ * D)) + #p = logistic(β₀ + β₁ * Sᵢ + β₂ * (S-Sᵢ) + β₃ * D) dist = Binomial(Sᵢ, p) - rand(dist) + draw = rand(dist) + + if return_prob + return draw, p + else + return draw + end end -function debris_transition(D, Y⃗, X⃗, δ=0.05, Γ=0.5) +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 @@ -42,16 +63,31 @@ function overall_losses(S⃗, D) #= Simulate the losses for each constellation =# - [ individual_loss(s, sum(S⃗),D) for s in S⃗ ] + + 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 @@ -75,70 +111,173 @@ Basic simulations 3. What happens when both participants maintain the cournot equilibrium? =# -#First try +#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 -policy_1(s,S,D) = 1 -function policy_2(s,S,D) - max(0,33-s) +struct SustainmentPolicy <: Policy + #= + Operator launches enough satellites to maintain a certain size + =# + level end -function policy_3(s,S,D) - max(0,20-s) +(p::SustainmentPolicy)(s,S,D) = max(0,p.level - s) + +struct simulation_result + state_path + profit_path + destruction_events + launch_choices + destruction_probabilities end -#setup state -n=200#00 -k=4 -state = fill([10.0,1,0,10],n) -profit_path = [] - -for i in 2:n - current_state = state[i-1] - #current_state = state[1] - S⃗ = current_state[2:end] - D = current_state[1] - - #Calculate policies - #X = [policy_2(s,sum(S⃗),D) for s in S⃗] - X = [2.0,2.0,policy_2(S⃗[3],0,0)] - #X = [policy_2(S⃗[1],sum(S⃗),D), policy_3(S⃗[2],sum(S⃗),D), ] - - #Information generation - price = inverse_demand_cournot(S⃗,10000) - push!(profit_path,[profit_function(s,x,1,2,0, 0,price) for (s,x) in zip(S⃗,X) ]) - - #calculate transition to next state - #loss of satellites - Y⃗ = overall_losses(S⃗, D) - - #debris - new_debris = debris_transition(D, Y⃗, X) - - #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) + +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) -debris = [v[1] for v in state] -stocks1 = [v[2] for v in state] -stocks2 = [v[3] for v in state] -stocks3 = [v[4] for v in state] plot( - [debris, stocks1, stocks2, stocks3] + simulation.state_path[2:4] + ,labels=["stocks 1" "stocks 2" "Stocks 3"] +) +plot( + simulation.state_path[1] + ,labels="debris" +) +plot( simulation.state_path ,labels=["debris" "stocks 1" "stocks 2" "Stocks 3"] ) -profit_1 = [v[1] for v in profit_path] -profit_2 = [v[2] for v in profit_path] -profit_3 = [v[3] for v in profit_path] +plot(simulation.profit_path) + +bar( + simulation.destruction_events + ,layout = grid(3,1) + ,legend=:none +) plot( - [profit_1,profit_2,profit_3] - ) \ No newline at end of file + 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 \ No newline at end of file