recording some past changes

CurrentWriting/Main.tex

@@ -53,14 +53,14 @@ and remaining policy questions.
     \subfile{sections/02_KesslerSyndrome}  %roughly done before 2021-07-14
     \subfile{sections/06_KesslerRegion}  %roughly done before 2021-07-14
 
-    \subsection{Constellation Operator's Program}\label{SEC:Operator}
+%     \subsection{Constellation Operator's Program}\label{SEC:Operator}
-    \subfile{sections/04_ConstellationOperator} %Reasonably done.
+%     \subfile{sections/04_ConstellationOperator} %Reasonably done.
-
+%
-    \subsection{Social Planner's Program}\label{SEC:Planner}
+%     \subsection{Social Planner's Program}\label{SEC:Planner}
-    \subfile{sections/05_SocialPlanner} %Reasonably done?
+%     \subfile{sections/05_SocialPlanner} %Reasonably done?
-
+%
-\section{Computation}\label{SEC:Computation}
+% \section{Computation}\label{SEC:Computation}
-\subfile{sections/07_ComputationalApproach} %needs some clarifications.
+% \subfile{sections/07_ComputationalApproach} %needs some clarifications.
 
 \section{Results}\label{SEC:Results}
 \subfile{sections/09_Results} %TODO

DecisionProcess.drawio

@@ -0,0 +1,87 @@
+<mxfile host="Electron" modified="2023-08-10T01:16:09.926Z" agent="Mozilla/5.0 (X11; Linux x86_64) AppleWebKit/537.36 (KHTML, like Gecko) draw.io/21.6.5 Chrome/114.0.5735.243 Electron/25.3.1 Safari/537.36" etag="ySVSO4vFxIwNG_GniQjf" version="21.6.5" type="device">
+  <diagram name="Page-1" id="1w1gsA99hJ9hG2iv5ZOv">
+    <mxGraphModel dx="1114" dy="843" grid="1" gridSize="10" guides="1" tooltips="1" connect="1" arrows="1" fold="1" page="1" pageScale="1" pageWidth="850" pageHeight="1100" math="0" shadow="0">
+      <root>
+        <mxCell id="0" />
+        <mxCell id="1" parent="0" />
+        <mxCell id="o47rVKH1999IgGybFU0q-11" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="o47rVKH1999IgGybFU0q-7" target="o47rVKH1999IgGybFU0q-10">
+          <mxGeometry relative="1" as="geometry" />
+        </mxCell>
+        <mxCell id="o47rVKH1999IgGybFU0q-21" style="rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=1;exitDx=0;exitDy=0;entryX=0.5;entryY=0;entryDx=0;entryDy=0;" edge="1" parent="1" source="o47rVKH1999IgGybFU0q-7" target="o47rVKH1999IgGybFU0q-15">
+          <mxGeometry relative="1" as="geometry" />
+        </mxCell>
+        <mxCell id="o47rVKH1999IgGybFU0q-34" style="rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=1;exitDx=0;exitDy=0;entryX=-0.012;entryY=0.5;entryDx=0;entryDy=0;entryPerimeter=0;" edge="1" parent="1" source="o47rVKH1999IgGybFU0q-7" target="o47rVKH1999IgGybFU0q-33">
+          <mxGeometry relative="1" as="geometry" />
+        </mxCell>
+        <mxCell id="o47rVKH1999IgGybFU0q-7" value="Environment State&lt;br&gt;&lt;ul&gt;&lt;li&gt;Each Constellation&#39;s&amp;nbsp; Size&lt;/li&gt;&lt;li&gt;Debris levels&lt;/li&gt;&lt;/ul&gt;" style="shape=process;whiteSpace=wrap;html=1;backgroundOutline=1;align=left;" vertex="1" parent="1">
+          <mxGeometry x="80" y="310" width="230" height="80" as="geometry" />
+        </mxCell>
+        <mxCell id="o47rVKH1999IgGybFU0q-12" style="rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;entryX=0;entryY=0;entryDx=0;entryDy=0;" edge="1" parent="1" source="o47rVKH1999IgGybFU0q-10" target="o47rVKH1999IgGybFU0q-9">
+          <mxGeometry relative="1" as="geometry" />
+        </mxCell>
+        <mxCell id="o47rVKH1999IgGybFU0q-22" style="rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=0.5;exitY=1;exitDx=0;exitDy=0;entryX=1;entryY=0;entryDx=0;entryDy=0;" edge="1" parent="1" source="o47rVKH1999IgGybFU0q-13" target="o47rVKH1999IgGybFU0q-15">
+          <mxGeometry relative="1" as="geometry" />
+        </mxCell>
+        <mxCell id="o47rVKH1999IgGybFU0q-23" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;" edge="1" parent="1" source="o47rVKH1999IgGybFU0q-15" target="o47rVKH1999IgGybFU0q-24">
+          <mxGeometry relative="1" as="geometry">
+            <mxPoint x="140" y="440" as="targetPoint" />
+          </mxGeometry>
+        </mxCell>
+        <mxCell id="o47rVKH1999IgGybFU0q-15" value="Environmental Transition Function" style="shape=dataStorage;whiteSpace=wrap;html=1;fixedSize=1;" vertex="1" parent="1">
+          <mxGeometry x="270" y="740" width="170" height="60" as="geometry" />
+        </mxCell>
+        <mxCell id="o47rVKH1999IgGybFU0q-25" style="rounded=0;orthogonalLoop=1;jettySize=auto;html=1;" edge="1" parent="1" source="o47rVKH1999IgGybFU0q-24" target="o47rVKH1999IgGybFU0q-7">
+          <mxGeometry relative="1" as="geometry" />
+        </mxCell>
+        <mxCell id="o47rVKH1999IgGybFU0q-24" value="Next Time Period" style="rhombus;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+          <mxGeometry x="155" y="470" width="80" height="80" as="geometry" />
+        </mxCell>
+        <mxCell id="o47rVKH1999IgGybFU0q-38" value="" style="group" vertex="1" connectable="0" parent="1">
+          <mxGeometry x="530" y="250" width="420" height="600" as="geometry" />
+        </mxCell>
+        <mxCell id="o47rVKH1999IgGybFU0q-32" value="Operators&#39; Problem" style="swimlane;whiteSpace=wrap;html=1;" vertex="1" parent="o47rVKH1999IgGybFU0q-38">
+          <mxGeometry x="20.005" y="9.677419354838714" width="367.5" height="600.0000000000001" as="geometry" />
+        </mxCell>
+        <mxCell id="o47rVKH1999IgGybFU0q-16" value="Profit" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="o47rVKH1999IgGybFU0q-32">
+          <mxGeometry x="32.49" y="490.32" width="157.5" height="59.68" as="geometry" />
+        </mxCell>
+        <mxCell id="o47rVKH1999IgGybFU0q-40" style="rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=0;exitY=0.25;exitDx=0;exitDy=0;entryX=1;entryY=0.5;entryDx=0;entryDy=0;" edge="1" parent="o47rVKH1999IgGybFU0q-32" source="o47rVKH1999IgGybFU0q-13" target="o47rVKH1999IgGybFU0q-16">
+          <mxGeometry relative="1" as="geometry" />
+        </mxCell>
+        <mxCell id="o47rVKH1999IgGybFU0q-13" value="Decision:&amp;nbsp;&lt;br&gt;&lt;span style=&quot;background-color: initial;&quot;&gt;How many satellites to launch&lt;/span&gt;" style="shape=hexagon;perimeter=hexagonPerimeter2;whiteSpace=wrap;html=1;fixedSize=1;align=center;" vertex="1" parent="o47rVKH1999IgGybFU0q-32">
+          <mxGeometry x="210.00166666666667" y="300.31967741935495" width="157.5" height="96.7741935483871" as="geometry" />
+        </mxCell>
+        <mxCell id="o47rVKH1999IgGybFU0q-9" value="Operator&#39;s State&lt;br&gt;&lt;ul&gt;&lt;li&gt;Constellation Size&lt;/li&gt;&lt;li&gt;Total Competitors&#39; Mass&lt;/li&gt;&lt;li&gt;Debris&lt;/li&gt;&lt;/ul&gt;" style="shape=process;whiteSpace=wrap;html=1;backgroundOutline=1;align=left;" vertex="1" parent="o47rVKH1999IgGybFU0q-32">
+          <mxGeometry x="60.00000000000001" y="90.32723180795199" width="226.04166666666669" height="153.2651162790698" as="geometry" />
+        </mxCell>
+        <mxCell id="o47rVKH1999IgGybFU0q-14" style="rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=0.5;exitY=1;exitDx=0;exitDy=0;" edge="1" parent="o47rVKH1999IgGybFU0q-32" source="o47rVKH1999IgGybFU0q-9" target="o47rVKH1999IgGybFU0q-13">
+          <mxGeometry relative="1" as="geometry" />
+        </mxCell>
+        <mxCell id="o47rVKH1999IgGybFU0q-41" style="rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=0.25;exitY=1;exitDx=0;exitDy=0;" edge="1" parent="o47rVKH1999IgGybFU0q-32" source="o47rVKH1999IgGybFU0q-9" target="o47rVKH1999IgGybFU0q-16">
+          <mxGeometry relative="1" as="geometry" />
+        </mxCell>
+        <mxCell id="o47rVKH1999IgGybFU0q-39" style="rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=1;exitDx=0;exitDy=0;exitPerimeter=0;entryX=0;entryY=0.5;entryDx=0;entryDy=0;" edge="1" parent="1" source="o47rVKH1999IgGybFU0q-33" target="o47rVKH1999IgGybFU0q-16">
+          <mxGeometry relative="1" as="geometry" />
+        </mxCell>
+        <mxCell id="o47rVKH1999IgGybFU0q-43" value="" style="group" vertex="1" connectable="0" parent="1">
+          <mxGeometry x="370" y="240" width="170" height="350" as="geometry" />
+        </mxCell>
+        <mxCell id="o47rVKH1999IgGybFU0q-10" value="Filter" style="triangle;whiteSpace=wrap;html=1;" vertex="1" parent="o47rVKH1999IgGybFU0q-43">
+          <mxGeometry x="80" y="70" width="60" height="80" as="geometry" />
+        </mxCell>
+        <mxCell id="o47rVKH1999IgGybFU0q-33" value="Market Mechanism&lt;br&gt;sets&amp;nbsp; price" style="shape=card;whiteSpace=wrap;html=1;" vertex="1" parent="o47rVKH1999IgGybFU0q-43">
+          <mxGeometry x="40" y="190" width="80" height="100" as="geometry" />
+        </mxCell>
+        <mxCell id="o47rVKH1999IgGybFU0q-42" value="Information Process" style="swimlane;whiteSpace=wrap;html=1;" vertex="1" parent="o47rVKH1999IgGybFU0q-43">
+          <mxGeometry width="170" height="350" as="geometry" />
+        </mxCell>
+        <mxCell id="o47rVKH1999IgGybFU0q-45" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;" edge="1" parent="1" source="o47rVKH1999IgGybFU0q-44" target="o47rVKH1999IgGybFU0q-15">
+          <mxGeometry relative="1" as="geometry" />
+        </mxCell>
+        <mxCell id="o47rVKH1999IgGybFU0q-44" value="Noise in Debris Generation" style="shape=document;whiteSpace=wrap;html=1;boundedLbl=1;" vertex="1" parent="1">
+          <mxGeometry x="295" y="890" width="120" height="80" as="geometry" />
+        </mxCell>
+      </root>
+    </mxGraphModel>
+  </diagram>
+</mxfile>

Notes_on_Model.txt

@@ -0,0 +1,70 @@
+# Thoughts on modelling
+
+Initial plan is to handle just two operators.
+
+## Environment
+
+### Transition function
+
+#### Debris
+- Currently using a decay formulation for debris.
+- Probably better to use a mass conversion from satellite -> debris
+    - This would have to be a stochastic approach to handle satellite destruction.
+    - I envision it would be better because it would represent how much debris 
+    is created better and would easigly lend itself to a size based debris 
+    approach.
+    - It doesn't represent slow degradation well though. It does handle destructive collisions better.
+- Maybe I can ignore satellite on satellite collisions in an initial model because it is the lowest.
+    
+#### Satellite Stocks
+
+- Probably just fine as is.
+    
+
+### Market
+
+- Current plan is to use a cournot market approach with inverse demand: \[ P = a-bQ \]
+- Need to look up dynamic cournot market lit.
+
+## Operator's Problem
+
+- Restrict policy to integers (e.g. 1:10)
+- Use a symmetric policy and value function for each operator.
+
+### Operator State
+- Initally I am going to track one type of debris, the operator's fleet size, and the number of satellites in all other fleets.
+
+## Social Planner's Problem
+
+I see two different classes of problems that the Social Planner could face, and both could take various forms.
+Here are the two areas and some variations:
+
+1. Maximizing long term expected utility/value
+    a. Maximize expected value
+    b. Preserve value x% of the time
+2. Preserving the resource (manage debris levels)
+    a. Keep debris below a threshold.
+    b. Keep debris from growing faster than Y rate.
+    c. Keep kessler syndrome from occuring in Z% over ZZ periods.
+    
+The first case "should" subsume the second in that minimal debris over a long period should permit more earnings. 
+But, if the time-discounting factor is low enough, this might not be the case.
+And this brings up a big issue in that the discount factor of the Social Planner and the Operators might differ 
+a lot.
+The question it poses is "Should we use this resource as we best see fit or preserve it for our descendants?"
+
+## modelling state
+
+I need to track all constellation sizes and debris size.
+
+# Thoughts on Interpretation etc,
+
+One option would be to take derivatives over the value function to get 
+Benefit functions. Look at when it is positive etc.
+
+Another thing to do would be to construct debris tracks.
+
+Another would be to search for steady states.
+
+I could compare different versions of the social planner's problem to see
+how they differ, i.e. in debris dynamics etc.

julia_code/attempt2/SimulationSteps.jl

@@ -0,0 +1,291 @@
+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

oci_container/Dockerfile

@@ -0,0 +1,25 @@
+FROM debian:latest
+MAINTAINER youainti@youainti.com
+RUN apt update
+RUN apt install -y vim
+RUN apt install -y curl
+
+#add user
+RUN useradd julia
+USER julia
+WORKDIR /home/julia
+#setup basic files
+RUN touch .profile
+RUN touch .bashrc
+
+#install julia
+RUN curl -fsSL https://install.julialang.org > juliaup_installer.sh
+RUN sh juliaup_installer.sh -y 
+RUN rm juliaup_installer.sh
+
+
+#install pluto
+RUN ~/.juliaup/bin/julia -e 'using Pkg; Pkg.add("Pluto")'
+
+#entrypoint
+CMD ~/.juliaup/bin/julia -e 'using Pluto; Pluto.run(;host="0.0.0.0")'