Optimization via Adaptive Resolvent Splitting
Date:
I described a novel method for building splitting algorithms for convex optimization problems. The method is based on the proximal resolvent operator, and it allows for the construction of adaptive splitting algorithms that can be tuned to the problem at hand.
Supporting notebook
Weapon Target Allocation Notebook
Use splitting to solve the weapon target allocation problem.
Author: Peter Barkley
Date: 12/14/2020
$$\min \sum_{i=1}^{n} V_i\prod_{j=1}^{m} (1 - Pk_{ij})^{x_{ij}}$$
$$\text{s.t.} \sum_{i=1}^{n} x_{ij} \le 1, \forall j$$
In [ ]:
# Libraries
import numpy as np
import cvxpy as cp
import pydeck as pdk
import pandas as pd
import wta
import geopy.distance
np.set_printoptions(precision=3, suppress=True)
In [ ]:
# Generate scenario
# five enemy combatants in an X,
# 6 friendly UxS: 2 UUV, 2 USV, 2 UAV
# Inventory: 6 torpedoes per UUV, 4 missiles per USV, 2 missiles per UAV
# Torpedoes: 2 nm range, high damage
# Missiles: 60 nm range from USV, medium damage
# Missiles: 80 nm range from UAV, medium damage
# Center enemy is highest value
# Front two enemy are lowest value
# Left back is medium value
# Right back is medium-high value
In [ ]:
# Location center
reference_location = np.array([23.5833, 119.5833])
ship_center = np.array([23.5833, 119.433])
uas_center = np.array([23.59, 119.9533])
In [ ]:
# Create json file for map with lat/long
# like [{"lat":50.5112014,"lon":6.9939896,"value":5,"type":"ferry"},
# Center enemy will be reference location plus 0.1 lat, 0.1 long
# Create json file for map with lat/long
n = 5
m = 6
# Lat and long each start with ship center and add offset to form X
tgt_lat = ship_center[0] + [0, 0.1, 0.1, -0.1, -0.1]
tgt_lon = ship_center[1] + [0, 0.1, -0.1, 0.1, -0.1]
# Lat and long each start with uas center and add offset to form X
uxs_lat = uas_center[0] + [0, 0.1, 0.1, -0.1, -0.1, -0.2]
uxs_lon = uas_center[1] + [0, 0.1, -0.1, 0.1, -0.1, 0]
lat = np.concatenate((tgt_lat, uxs_lat))
lon = np.concatenate((tgt_lon, uxs_lon))
tgt_values = 10*np.array([5, 1, 1, 3, 4])
value = np.concatenate((tgt_values, [0, 0, 0, 0, 0, 0]))
shiptypes = np.array(["ferry", "ship", "ship", "ferry", "ferry"])
uxvtypes = np.array(["uuv", "uuv", "usv", "usv", "uav", "uav"])
platform_types = np.concatenate((shiptypes, uxvtypes))
# data = np.concatenate(np.array([lat, lon, value]), shiptypes)
# data = data.T
# data = data.tolist()
# data = [{"lat": data[i][0], "lon": data[i][1], "value": data[i][2], "type": data[i][3]} for i in range(n+m)]
# # Save data to json file
# import json
# with open("data.json", "w") as outfile:
# json.dump(data, outfile)
# # Print data
# print(data)
In [ ]:
# Map
BASE_URL = 'https://raw.githubusercontent.com/peterbarkley/wta/main/images/'
FERRY_URL = BASE_URL + 'ferry_right.png'
SHIP_URL = BASE_URL + 'warship_054A_down_right.png'
SUB_URL = BASE_URL + 'sub_up_left.png'
UAV_URL = BASE_URL + 'drone_icon_up_left.png'
USV_URL = BASE_URL + 'ship_icon_up_left.png'
# Dictionary of URLs by platform
ICON_URL = {
"ferry": FERRY_URL,
"ship": SHIP_URL,
"sub": SUB_URL,
"uav": UAV_URL,
"uuv": SUB_URL,
"usv": USV_URL,
}
dim = {}
for platform in ICON_URL:
if platform in ["uav", "uuv", "usv"]:
dim[platform] = 120
else:
dim[platform] = 242
# Icon data for each platform (UAV, FERRY, SHIP, SUB)
icon_data = {}
for platform in ICON_URL:
d = dim[platform]
icon_data[platform] = {
"url": ICON_URL[platform],
"width": d,
"height": d,
#"anchorY": 0,
}
data = pd.read_json("data.json")
data["icon_data"] = data["type"].apply(lambda x: icon_data[x])
In [ ]:
# Map data for values map
enemy_data = data[data["value"] > 0]
friendly_data = data[data["value"] == 0]
view_state = pdk.ViewState(
latitude=float(reference_location[0]),
longitude=float(reference_location[1]),
zoom=10,
bearing=-45,
pitch=60,
)
icon_layer = pdk.Layer(
type="IconLayer",
data=data,
get_icon="icon_data",
get_size=4,
size_scale=15,
get_position=["lon", "lat"],
pickable=True,
auto_highlight=True,
)
column_layer = pdk.Layer(
"ColumnLayer",
data=enemy_data,
get_position=["lon", "lat"],
get_elevation="value",
elevation_scale=500,
radius=500,
get_fill_color=["value * 100", "value", "value * 100", 140],
pickable=True,
auto_highlight=True,
)
tooltip = {
"html": "Value <b>{value}</b>",
"style": {"background": "grey", "color": "white", "font-family": '"Helvetica Neue", Arial', "z-index": "10000"},
}
r = pdk.Deck(layers=[icon_layer, column_layer],
tooltip=tooltip,
initial_view_state=view_state,
map_style=pdk.map_styles.LIGHT,)# map_provider="mapbox", map_style=pdk.map_styles.SATELLITE)
r.to_html("icon_layer.html")
Out[ ]:
In [ ]:
#Export as image
# import imgkit
# imgkit.from_file('icon_layer.html', 'toy_WTA_values.jpg')
In [ ]:
# Generate scenario data for target values, Pk, and weapon's per UxS
uxs_weapons = {'uav': 2, 'uuv': 4, 'usv': 6}
uxs_base_pk = {'uav': 0.2, 'uuv': 0.3, 'usv': 0.1}
weapon_count = [uxs_weapons[t] for t in uxvtypes]
q = np.zeros((n,m))
ds = []
for i in range(n):
for j in range(m):
d = geopy.distance.distance((tgt_lat[i], tgt_lon[i]), (uxs_lat[j], uxs_lon[j])).nm
loc_damage = 1/(d-2)
ds.append(d)
q[i, j] = 1 - (min(.05, loc_damage) + uxs_base_pk[uxvtypes[j]])+ 0.1/tgt_values[i]
print(q)
[[0.665 0.671 0.857 0.871 0.756 0.767] [0.664 0.672 0.86 0.875 0.76 0.773] [0.679 0.683 0.872 0.885 0.775 0.783] [0.658 0.669 0.853 0.866 0.753 0.758] [0.672 0.677 0.868 0.876 0.765 0.772]]
In [ ]:
# Build arc layer showing Pk
# from q ranging from almost white for the max q to full red for the min q
q_min = np.min(q)
q_max = np.max(q)
# Scale q to range from 0 to 1
q_scaled = (q - q_min)/(q_max - q_min)
# Scale q to range from 0 to 255
q_scaled = q_scaled * 255
# Convert q to integer
q_scaled = q_scaled.astype(int)
# Scale weapon count
weapon_count_scaled = np.array(weapon_count)
weapon_count_scaled = (np.max(weapon_count_scaled)-weapon_count_scaled)/np.max(weapon_count_scaled)
weapon_count_scaled = weapon_count_scaled * 255
weapon_count_scaled = weapon_count_scaled.astype(int)
#print(weapon_count_scaled)
# Create data for arc layer
arc_data = []
for i in range(n):
for j in range(m):
arc_data.append([lon[i], lat[i], lon[j+n], lat[j+n], q_scaled[i,j], round(1-q[i,j],2), weapon_count_scaled[j], weapon_count[j]])
# Create data for arc layer
arc_df = pd.DataFrame(arc_data, columns=["lon", "lat", "lon2", "lat2", "qs", "Pk", "wcs", "weapon_count"])
# Create arc layer
arc_layer = pdk.Layer(
type="ArcLayer",
data=arc_df,
get_source_position=["lon2", "lat2"],
get_target_position=["lon", "lat"],
get_source_color=[255, 255, 255], # want to use [255, "wcs", 255] but not working without manual switch
get_target_color=["qs", 0, 0],
get_width=2,
pickable=True,
auto_highlight=True,
)
tooltip = {
"html": "Pk <b>{Pk}</b>, Weapon Count <b>{weapon_count}</b>",
"style": {"background": "grey", "color": "white", "font-family": '"Helvetica Neue", Arial', "z-index": "10000"},
}
# Create deck
r = pdk.Deck(layers=[icon_layer, arc_layer],
tooltip=tooltip,
initial_view_state=view_state,
map_style=pdk.map_styles.LIGHT,)
# Save deck to html file
r.to_html("pk_arc_layer.html")
Out[ ]:
In [ ]:
# Function to get final probability of kill for each target
def get_final_surv_prob(q, x):
"""
Get the final probability of kill for each target.
Inputs:
q: (n,m) array of survival probabilities
x: (n,m) array of weapon assignments
"""
return np.prod(np.power(q, x), axis=1)
In [ ]:
# Function to get the value of the problem if each platform solves independently
def get_ind_value(q, V, W):
"""
Get the total value if each platform solves independently.
Inputs:
q: (n,m) array of survival probabilities
V: (n,) array of target values
W: (m,) array of weapon counts
"""
# Loop through platforms
n, m = q.shape
x = np.zeros((n,m))
for i in range(m):
# Solve the WTA problem for platform i
q_i = q[:,i]
q_i = q_i.reshape((n,1))
pv, x_i = wta.wta(q_i, V, W[i])
x[:,i] = x_i[:,0]
return V@get_final_surv_prob(q, x), x
In [ ]:
# Solve independent WTA problem
V = tgt_values
WW = weapon_count
ind_value, ind_x = get_ind_value(q, V, WW)
ind_surv_prob = get_final_surv_prob(q, ind_x)
print(ind_value)
print(ind_x)
print(ind_surv_prob)
optimal optimal optimal optimal optimal optimal 41.22014561081071 [[2. 2. 4. 4. 1. 1.] [0. 0. 0. 0. 0. 0.] [0. 0. 0. 0. 0. 0.] [1. 1. 0. 0. 0. 0.] [1. 1. 2. 2. 1. 1.]] [0.036 1. 1. 0.44 0.155]
In [ ]:
# Function to display icons and columns and arc based on arcs implied by solution
def getSolutionArcMap(x, surv_prob, q, q_scaled, enemy_data, friendly_data, weapon_count):
"""
Get the arcs implied by the solution.
Inputs:
x: (n,m) array of weapon assignments
enemy_data: (n,3) dataframe of enemy data
friendly_data: (m,3) dataframe of friendly data
weapon_count: (m,) array of weapon counts
"""
n, m = x.shape
# Create data for arc layer
arc_data = []
for i in range(n):
for j in range(m):
if x[i,j] > 0.5:
e_lon = enemy_data["lon"][i]
e_lat = enemy_data["lat"][i]
f_lon = friendly_data["lon"][n+j]
f_lat = friendly_data["lat"][n+j]
arc_data.append([e_lon, e_lat, f_lon, f_lat, q_scaled[i,j], round(1-q[i,j],2), round(x[i,j]), weapon_count[j]])
# Create data for arc layer
arc_df = pd.DataFrame(arc_data, columns=["lon", "lat", "lon2", "lat2", "qs", "Pk", "used", "avail"])
# Create arc layer
arc_layer = pdk.Layer(
type="ArcLayer",
data=arc_df,
get_source_position=["lon2", "lat2"],
get_target_position=["lon", "lat"],
get_source_color=[255, 255, 255], # want to use [255, "wcs", 255] but not working without manual switch
get_target_color=["qs", 0, 0],
get_width=2,
pickable=True,
auto_highlight=True,
)
tooltip = {
"html": "Pk <b>{Pk}</b>, Used <b>{used}</b>, Available <b>{avail}</b>",
"style": {"background": "grey", "color": "white", "font-family": '"Helvetica Neue", Arial', "z-index": "10000"},
}
# Create data for column layer
column_data = []
for i in range(n):
e_lon = enemy_data["lon"][i]
e_lat = enemy_data["lat"][i]
tgt_val = enemy_data["value"][i]
column_data.append([e_lon, e_lat, tgt_val*surv_prob[i]])
# Create data for column layer
column_df = pd.DataFrame(column_data, columns=["lon", "lat", "value"])
# Create column layer
column_layer = pdk.Layer(
"ColumnLayer",
data=column_df,
get_position=["lon", "lat"],
get_elevation="value",
elevation_scale=500,
radius=500,
get_fill_color=["value * 100", "value", "value * 100", 140],
pickable=True,
auto_highlight=True,
)
# Tooltip
tooltip_col = {
"html": "Value <b>{value}</b>",
"style": {"background": "grey", "color": "white", "font-family": '"Helvetica Neue", Arial', "z-index": "10000"},
}
# Create deck
r = pdk.Deck(layers=[icon_layer, arc_layer, column_layer],
tooltip=tooltip,
initial_view_state=view_state,
map_style=pdk.map_styles.LIGHT,)
# Save deck to html file
return r
In [ ]:
# Build independent arc layer map
r = getSolutionArcMap(ind_x, ind_surv_prob, q, q_scaled, enemy_data, friendly_data, weapon_count)
r.to_html("ind_arc_layer.html")
Out[ ]:
In [ ]:
# Solve full WTA problem for both integer and continuous values
prob_val_int, x_full_int = wta.wta(q, tgt_values, weapon_count)
print("Integer")
print(prob_val_int)
print("Non-integer")
prov_val, x_full = wta.wta(q, tgt_values, weapon_count, integer=False)
print(prov_val)
optimal Integer 33.68999939376743 Non-integer optimal 33.628945542954824
In [ ]:
# Get map for integer solution
r = getSolutionArcMap(x_full_int, get_final_surv_prob(q, x_full_int), q, q_scaled, enemy_data, friendly_data, weapon_count)
r.to_html("integer_arc_layer.html")
Out[ ]:
In [ ]:
# Get map for non-integer solution
r = getSolutionArcMap(x_full, get_final_surv_prob(q, x_full), q, q_scaled, enemy_data, friendly_data, weapon_count)
r.to_html("non_integer_arc_layer.html")
Out[ ]:
In [ ]:
# Get map with 7000 enemy platforms and 3000 friendly platforms
# Spread to the north and south of the reference location
# Copy enemy data
huge_enemy_data = enemy_data.copy()
# Loop over data shifting lat and lon
for i in range(7000):
huge_enemy_data.loc[i] = enemy_data.loc[i%n]
huge_enemy_data.loc[i, "lat"] = reference_location[0] - 1.1 + i/1500 + np.random.normal(0, 0.2)
huge_enemy_data.loc[i, "lon"] = reference_location[1] -1/60 + i/5000 + min(np.random.normal(0, 0.15), .15)
In [ ]:
# Copy friendly data
huge_friendly_data = friendly_data.copy()
# Loop over data shifting lat and lon
for i in range(200):
huge_friendly_data.loc[i] = friendly_data.loc[n + i%m]
huge_friendly_data.loc[i, "lat"] = reference_location[0] - 2 + i/50 + np.random.normal(0, 0.25)
huge_friendly_data.loc[i, "lon"] = reference_location[1] + 3 + i/150 + np.random.normal(0, 0.2)
# Create icon layer
icon_layer_enemy_h = pdk.Layer(
type="IconLayer",
data=huge_enemy_data,
get_icon="icon_data",
get_size=4,
size_scale=15,
get_position=["lon", "lat"],
pickable=True,
auto_highlight=True,
)
icon_layer_friendly_h = pdk.Layer(
type="IconLayer",
data=huge_friendly_data,
get_icon="icon_data",
get_size=4,
size_scale=15,
get_position=["lon", "lat"],
pickable=True,
auto_highlight=True,
)
# # Create column layer
# column_layer = pdk.Layer(
# "ColumnLayer",
# data=huge_enemy_data,
# get_position=["lon", "lat"],
# get_elevation="value",
# elevation_scale=500,
# radius=500,
# get_fill_color=["value * 100", "value", "value * 100", 140],
# pickable=True,
# auto_highlight=True,
# )
In [ ]:
# Generate map for huge scenario
# Tooltip
tooltip = {
"html": "Value <b>{value}</b>",
"style": {"background": "grey", "color": "white", "font-family": '"Helvetica Neue", Arial', "z-index": "10000"},
}
# Create map
r = pdk.Deck(layers=[icon_layer_enemy_h, icon_layer_friendly_h],
tooltip=tooltip,
initial_view_state=view_state,
map_style=pdk.map_styles.LIGHT,)
# Save deck to html file
r.to_html("huge_icon_layer.html")
Out[ ]:
In [ ]:
# Get map with 120 enemy platforms and 20 friendly platforms
# Spread to the north and south of the reference location
medium_n = 120
medium_m = 20
# Copy enemy data
medium_enemy_data = enemy_data.copy()
# Loop over data shifting lat and lon
for i in range(medium_n):
medium_enemy_data.loc[i] = enemy_data.loc[i%n]
medium_enemy_data.loc[i, "lat"] = reference_location[0] - 1.7 + i/150 + np.random.normal(0, 0.2)
medium_enemy_data.loc[i, "lon"] = reference_location[1] - 1/60 + i/500 + min(np.random.normal(0, 0.15), .15)
In [ ]:
# Copy friendly data
medium_friendly_data = friendly_data.copy()
# Loop over data shifting lat and lon
for i in range(medium_m):
medium_friendly_data.loc[i] = friendly_data.loc[n + i%m]
medium_friendly_data.loc[i, "lat"] = reference_location[0] - 1.5 + i/50 + np.random.normal(0, 0.1)
medium_friendly_data.loc[i, "lon"] = reference_location[1] + 0.5 + i/100 + np.random.normal(0, 0.1)
# Create icon layer
icon_layer_enemy_medium = pdk.Layer(
type="IconLayer",
data=medium_enemy_data,
get_icon="icon_data",
get_size=4,
size_scale=15,
get_position=["lon", "lat"],
pickable=True,
auto_highlight=True,
)
icon_layer_friendly_medium = pdk.Layer(
type="IconLayer",
data=medium_friendly_data,
get_icon="icon_data",
get_size=4,
size_scale=15,
get_position=["lon", "lat"],
pickable=True,
auto_highlight=True,
)
# # Create column layer
# column_layer = pdk.Layer(
# "ColumnLayer",
# data=huge_enemy_data,
# get_position=["lon", "lat"],
# get_elevation="value",
# elevation_scale=500,
# radius=500,
# get_fill_color=["value * 100", "value", "value * 100", 140],
# pickable=True,
# auto_highlight=True,
# )
In [ ]:
# Import solution from hamming (alg_x_120_20.json)
import json
with open("alg_x_120_20.json", "r") as infile:
alg_x = json.load(infile)
In [ ]:
# Create arc layer for medium scenario
medium_arc_data = []
for i in range(medium_n):
for j in range(medium_m):
if alg_x[i][j] > 0.5:
e_lon = medium_enemy_data["lon"][i]
e_lat = medium_enemy_data["lat"][i]
f_lon = medium_friendly_data["lon"][j]
f_lat = medium_friendly_data["lat"][j]
medium_arc_data.append([e_lon, e_lat, f_lon, f_lat])
# Create data for arc layer
medium_arc_df = pd.DataFrame(medium_arc_data, columns=["lon", "lat", "lon2", "lat2"])
# Create arc layer
arc_layer_medium = pdk.Layer(
type="ArcLayer",
data=medium_arc_df,
get_source_position=["lon2", "lat2"],
get_target_position=["lon", "lat"],
get_source_color=[255, 255, 255], # want to use [255, "wcs", 255] but not working without manual switch
get_target_color=[255, 0, 0],
get_width=2,
pickable=False,
auto_highlight=True,
)
In [ ]:
# Get final values for medium scenario solution
med_q, med_V, med_W = wta.generate_random_problem(120, 20)
med_surv_prob = get_final_surv_prob(med_q, alg_x)
med_value = np.multiply(med_V, med_surv_prob)
In [ ]:
# Generate column data for medium scenario
column_data = []
for i in range(medium_n):
e_lon = medium_enemy_data["lon"][i]
e_lat = medium_enemy_data["lat"][i]
tgt_val = med_value[i]
column_data.append([e_lon, e_lat, tgt_val])
# Create data for column layer
column_df = pd.DataFrame(column_data, columns=["lon", "lat", "value"])
# Create column layer
column_layer_medium = pdk.Layer(
"ColumnLayer",
data=column_df,
get_position=["lon", "lat"],
get_elevation="value",
elevation_scale=500,
radius=500,
get_fill_color=["value * 100", "value", "value * 100", 140],
pickable=False,
auto_highlight=True,
)
In [ ]:
# Build medium scenario map
# Tooltip
tooltip = {
"html": "Value <b>{value}</b>",
"style": {"background": "grey", "color": "white", "font-family": '"Helvetica Neue", Arial', "z-index": "10000"},
}
# Create map
r = pdk.Deck(layers=[icon_layer_enemy_medium,
icon_layer_friendly_medium,
arc_layer_medium,
column_layer_medium],
initial_view_state=view_state,
map_style=pdk.map_styles.LIGHT,)
# Save deck to html file
r.to_html("medium_icon_layer.html")
Out[ ]:
In [ ]:
# Get map for permissible communication paths as arcs
# Platform 6 can only talk to platforms 9 and 10
# Create data for arc layer
comms_arc_data = []
for j in range(m):
f_lon = friendly_data["lon"][n+j]
f_lat = friendly_data["lat"][n+j]
if n+j == 6:
for k in [9, 10]:
f_lon2 = friendly_data["lon"][k]
f_lat2 = friendly_data["lat"][k]
comms_arc_data.append([f_lon, f_lat, f_lon2, f_lat2])
else:
for k in range(j+1, m):
if n+k != 6:
f_lon2 = friendly_data["lon"][n+k]
f_lat2 = friendly_data["lat"][n+k]
comms_arc_data.append([f_lon, f_lat, f_lon2, f_lat2])
In [ ]:
# Generate map for permissible communication paths
# Create data for arc layer
comms_arc_df = pd.DataFrame(comms_arc_data, columns=["lon", "lat", "lon2", "lat2"])
# Create arc layer
comms_arc_layer = pdk.Layer(
type="ArcLayer",
data=comms_arc_df,
get_source_position=["lon2", "lat2"],
get_target_position=["lon", "lat"],
get_source_color=[255, 255, 255],
get_target_color=[255, 255, 255],
get_width=2,
pickable=True,
auto_highlight=True,
)
r = pdk.Deck(layers=[icon_layer, comms_arc_layer],
initial_view_state=view_state,
map_style=pdk.map_styles.LIGHT,)
r.to_html("comm_arc_layer.html")
Out[ ]:
In [ ]:
# Algorithm matrices
M = np.array([[-1, 1, 0, 0, 0, 0],
[0, -1, 1, 0, 0, 0],
[0, 0, -1, 1, 0, 0],
[0, 0, 0, -1, 1, 0],
[0, 0, 0, 0, -1, 1]])
W = M.T@M
L = np.array([[0, 0, 0, 0, 0, 0],
[1, 0, 0, 0, 0, 0],
[0, 1, 0, 0, 0, 0],
[0, 0, 1, 0, 0, 0],
[0, 0, 0, 1, 0, 0],
[1, 0, 0, 0, 1, 0]])
# print(np.linalg.eigvals(W + L + L.T - 2*np.eye(6)))
print(W)
print(L)
[[ 1 -1 0 0 0 0] [-1 2 -1 0 0 0] [ 0 -1 2 -1 0 0] [ 0 0 -1 2 -1 0] [ 0 0 0 -1 2 -1] [ 0 0 0 0 -1 1]] [[0 0 0 0 0 0] [1 0 0 0 0 0] [0 1 0 0 0 0] [0 0 1 0 0 0] [0 0 0 1 0 0] [1 0 0 0 1 0]]
In [ ]:
# Data
n = 6
tgts = 5
wpns = 6
m = (tgts, wpns)
# Target values
V = tgt_values
# Weapon values
WW = weapon_count
In [ ]:
# Verify solution
p, x = wta.wta(q, V, WW, integer=False)
optimal
In [ ]:
# Algorithm initialization
itrs = 501
gamma = 0.5
v0 = []
vk = []
log = []
log_e = []
for i in range(n):
vi = 1/tgts*np.array(WW)*np.ones(m)
v0.append(vi)
vk.append(vi.copy())
node_tgts = {0:[0, 1], 1:[1,2], 2:[2,3], 3:[3,4], 4:[4,0], 5:[1, 4]}
num_nodes_per_tgt = [2, 3, 2, 2, 3]
In [ ]:
# Create variables/params/objs for the algorithm
probs = [] # List of problems for each node
all_x = [] # List of last x solutions for each node as params
all_v = [] # List of last v solutions for each node as params
all_w = []# List of variables for each node
for i in range(n):
w = cp.Variable(m)
all_w.append(w)
x = cp.Parameter(m)
x.value = np.zeros(m)
all_x.append(x)
v = cp.Parameter(m)
v.value = v0[i]
all_v.append(v)
y = v + sum(L[i,j]*all_x[j] for j in range(i)) + L[i,i]*w
qq = np.ones(m)
qq[node_tgts[i]] = q[node_tgts[i]] # Only use the targets that are in the node
weighted_weapons = cp.multiply(w, np.log(qq)) # (tgts, wpns)
survival_probs = cp.exp(cp.sum(weighted_weapons, axis=1)) # (tgts,)
VV = np.zeros(tgts)
VV[node_tgts[i]] = V[node_tgts[i]]
VV = VV/num_nodes_per_tgt
obj = cp.Minimize(VV@survival_probs + .5*cp.sum_squares(w - y))
cons = [w >= 0, cp.sum(w, axis=0) <= WW]
probs.append(cp.Problem(obj, cons))
In [ ]:
# Run algorithm
for itr in range(itrs):
#print("Iteration", itr)
e = 0
for i in range(n):
probs[i].solve()
if itr % 100 == 0:
print("Iteration", itr, "Node", i)
print(all_w[i].value)
e += np.linalg.norm(all_w[i].value - all_x[i].value)
all_x[i].value = all_w[i].value
log_e.append(e)
for i in range(n):
vk[i] -= gamma*sum(W[i,j]*all_x[j].value for j in range(n))
all_v[i].value = vk[i]
log.append(V@get_final_surv_prob(q, all_x[0].value))
if itr % 100 == 0:
print("v", vk)
Iteration 0 Node 0
[[1.603 1.583 1.505 1.472 0.95 0.922]
[0.88 0.876 1.227 1.223 0.451 0.446]
[0.506 0.514 1.089 1.102 0.2 0.211]
[0.506 0.514 1.089 1.102 0.2 0.211]
[0.506 0.514 1.089 1.102 0.2 0.211]]
Iteration 0 Node 1
[[1.428 1.411 1.442 1.415 0.834 0.811]
[1.019 1.009 1.28 1.269 0.545 0.533]
[0.894 0.897 1.225 1.224 0.454 0.457]
[0.33 0.341 1.026 1.045 0.083 0.1 ]
[0.33 0.341 1.026 1.045 0.083 0.1 ]]
Iteration 0 Node 2
[[1.1 1.095 1.321 1.305 0.58 0.568]
[0.692 0.693 1.159 1.159 0.292 0.29 ]
[1.047 1.053 1.273 1.267 0.516 0.519]
[1.159 1.135 1.343 1.333 0.612 0.623]
[0.003 0.025 0.905 0.935 0. 0. ]]
Iteration 0 Node 3
[[0.684 0.689 1.169 1.165 0.299 0.295]
[0.275 0.288 1.007 1.019 0.011 0.017]
[0.63 0.648 1.122 1.126 0.235 0.246]
[1.493 1.45 1.476 1.451 0.839 0.848]
[0.918 0.925 1.226 1.239 0.616 0.594]]
Iteration 0 Node 4
[[1.51 1.498 1.501 1.457 0.779 0.755]
[0. 0. 0.848 0.875 0. 0. ]
[0.16 0.193 0.962 0.982 0. 0. ]
[1.023 0.996 1.317 1.307 0.432 0.47 ]
[1.307 1.313 1.372 1.38 0.789 0.775]]
Iteration 0 Node 5
[[1.914 1.901 1.7 1.632 1.012 0.995]
[0.366 0.361 1.022 1.025 0.193 0.195]
[0. 0. 0.746 0.787 0. 0. ]
[0.329 0.329 1.1 1.112 0. 0. ]
[1.391 1.409 1.432 1.445 0.795 0.81 ]]
v [array([[0.712, 0.714, 1.168, 1.172, 0.342, 0.344],
[0.87 , 0.867, 1.227, 1.223, 0.447, 0.443],
[0.994, 0.992, 1.268, 1.261, 0.527, 0.523],
[0.712, 0.714, 1.168, 1.172, 0.342, 0.344],
[0.712, 0.714, 1.168, 1.172, 0.342, 0.344]]), array([[0.724, 0.728, 1.171, 1.173, 0.331, 0.334],
[0.567, 0.575, 1.113, 1.122, 0.226, 0.235],
[0.682, 0.687, 1.156, 1.16 , 0.304, 0.308],
[1.302, 1.283, 1.39 , 1.372, 0.722, 0.717],
[0.724, 0.728, 1.171, 1.173, 0.417, 0.406]]), array([[0.755, 0.756, 1.185, 1.185, 0.386, 0.385],
[0.755, 0.756, 1.185, 1.185, 0.386, 0.385],
[0.515, 0.519, 1.1 , 1.109, 0.228, 0.233],
[0.553, 0.561, 1.108, 1.115, 0.249, 0.251],
[1.421, 1.408, 1.422, 1.407, 0.75 , 0.747]]), array([[1.421, 1.407, 1.442, 1.416, 0.781, 0.767],
[0.871, 0.859, 1.196, 1.198, 0.535, 0.528],
[0.773, 0.775, 1.196, 1.198, 0.423, 0.414],
[0.398, 0.415, 1.054, 1.069, 0.083, 0.098],
[0.537, 0.544, 1.112, 1.119, 0.178, 0.193]]), array([[0.589, 0.597, 1.134, 1.141, 0.276, 0.289],
[1.121, 1.124, 1.367, 1.347, 0.502, 0.506],
[0.955, 0.931, 1.171, 1.175, 0.518, 0.523],
[0.688, 0.694, 1.171, 1.175, 0.387, 0.355],
[0.648, 0.654, 1.157, 1.161, 0.317, 0.327]]), array([[0.598, 0.599, 1.1 , 1.112, 0.284, 0.28 ],
[0.617, 0.619, 1.113, 1.125, 0.303, 0.302],
[0.88 , 0.897, 1.308, 1.297, 0.4 , 0.4 ],
[1.147, 1.133, 1.308, 1.297, 0.616, 0.635],
[0.758, 0.752, 1.17 , 1.168, 0.397, 0.382]])]
Iteration 100 Node 0
[[1.367 1.391 1.752 1.553 0.917 0.58 ]
[0.198 0.118 1.029 0.839 0. 0. ]
[0.143 0.378 0.533 0.56 0. 0. ]
[1.006 0.634 1.618 1.628 0.482 0.775]
[1.286 1.479 1.068 1.42 0.601 0.645]]
Iteration 100 Node 1
[[1.366 1.39 1.758 1.556 0.919 0.576]
[0.199 0.117 1.03 0.836 0. 0. ]
[0.143 0.384 0.525 0.552 0. 0. ]
[1.006 0.626 1.624 1.635 0.482 0.78 ]
[1.286 1.483 1.062 1.421 0.599 0.643]]
Iteration 100 Node 2
[[1.365 1.389 1.764 1.559 0.922 0.572]
[0.201 0.116 1.032 0.833 0. 0. ]
[0.143 0.389 0.517 0.544 0. 0. ]
[1.006 0.619 1.631 1.642 0.481 0.786]
[1.286 1.487 1.056 1.422 0.597 0.642]]
Iteration 100 Node 3
[[1.363 1.389 1.77 1.562 0.925 0.568]
[0.202 0.115 1.033 0.83 0. 0. ]
[0.143 0.395 0.51 0.536 0. 0. ]
[1.005 0.611 1.638 1.649 0.48 0.791]
[1.286 1.49 1.049 1.423 0.595 0.641]]
Iteration 100 Node 4
[[1.362 1.388 1.776 1.564 0.927 0.564]
[0.203 0.114 1.035 0.827 0. 0. ]
[0.143 0.401 0.502 0.528 0. 0. ]
[1.005 0.603 1.644 1.656 0.48 0.796]
[1.287 1.494 1.043 1.424 0.593 0.64 ]]
Iteration 100 Node 5
[[1.361 1.387 1.782 1.567 0.93 0.56 ]
[0.205 0.113 1.036 0.824 0. 0. ]
[0.143 0.407 0.494 0.521 0. 0. ]
[1.005 0.595 1.651 1.663 0.479 0.801]
[1.287 1.498 1.037 1.425 0.591 0.639]]
v [array([[ 0.452, 0.499, 1.407, 1.245, 0.291, -0.018],
[-0.262, -0.328, 0.863, 0.691, -0.305, -0.284],
[ 0.602, 0.827, 0.701, 0.708, 0.311, 0.293],
[ 1.464, 1.076, 1.793, 1.784, 0.793, 1.071],
[ 1.744, 1.927, 1.237, 1.573, 0.911, 0.937]]), array([[ 1.258, 1.248, 1.371, 1.35 , 0.707, 0.686],
[ 0.341, 0.357, 1.027, 1.045, 0.088, 0.112],
[-0.116, -0.099, 0.872, 0.903, -0.202, -0.183],
[ 1.259, 1.241, 1.371, 1.354, 0.704, 0.695],
[ 1.259, 1.253, 1.358, 1.348, 0.703, 0.689]]), array([[ 1.349, 1.334, 1.407, 1.384, 0.77 , 0.753],
[ 1.351, 1.334, 1.403, 1.379, 0.767, 0.757],
[-0.025, -0.013, 0.908, 0.938, -0.138, -0.116],
[-0.025, 0.007, 0.887, 0.916, -0.164, -0.15 ],
[ 1.35 , 1.338, 1.395, 1.383, 0.765, 0.756]]), array([[ 1.257, 1.243, 1.375, 1.358, 0.712, 0.698],
[ 1.26 , 1.243, 1.371, 1.353, 0.71 , 0.702],
[ 1.258, 1.25 , 1.362, 1.348, 0.709, 0.702],
[-0.117, -0.084, 0.855, 0.89 , -0.221, -0.205],
[ 0.341, 0.348, 1.037, 1.051, 0.09 , 0.103]]), array([[-0.115, -0.09 , 0.856, 0.894, -0.224, -0.195],
[ 1.259, 1.246, 1.371, 1.351, 0.711, 0.697],
[ 1.258, 1.253, 1.361, 1.346, 0.711, 0.697],
[ 1.258, 1.239, 1.376, 1.36 , 0.711, 0.702],
[ 0.341, 0.352, 1.037, 1.049, 0.091, 0.099]]), array([[ 0.6 , 0.567, 0.784, 0.969, 0.144, 0.476],
[ 0.851, 0.948, 1.166, 1.381, 0.429, 0.417],
[ 1.824, 1.583, 1.996, 1.957, 1.009, 1.007],
[ 0.961, 1.32 , 0.918, 0.896, 0.577, 0.285],
[-0.236, -0.418, 1.136, 0.796, -0.16 , -0.185]])]
Iteration 200 Node 0
[[1.303 1.349 2.054 1.687 1.044 0.376]
[0.27 0.073 1.102 0.694 0. 0. ]
[0.134 0.661 0.144 0.171 0. 0. ]
[0.997 0.242 1.944 1.975 0.451 1.034]
[1.297 1.674 0.755 1.473 0.505 0.59 ]]
Iteration 200 Node 1
[[1.302 1.348 2.061 1.689 1.047 0.372]
[0.271 0.072 1.103 0.692 0. 0. ]
[0.134 0.667 0.136 0.163 0. 0. ]
[0.997 0.235 1.951 1.982 0.451 1.039]
[1.297 1.678 0.749 1.475 0.503 0.589]]
Iteration 200 Node 2
[[1.3 1.347 2.067 1.692 1.049 0.368]
[0.273 0.071 1.105 0.689 0. 0. ]
[0.134 0.673 0.129 0.155 0. 0. ]
[0.996 0.227 1.957 1.989 0.45 1.044]
[1.297 1.682 0.743 1.476 0.501 0.588]]
Iteration 200 Node 3
[[1.299 1.346 2.073 1.695 1.052 0.364]
[0.274 0.07 1.106 0.686 0. 0. ]
[0.133 0.679 0.121 0.147 0. 0. ]
[0.996 0.219 1.964 1.996 0.449 1.049]
[1.297 1.686 0.736 1.477 0.499 0.587]]
Iteration 200 Node 4
[[1.298 1.345 2.079 1.697 1.054 0.36 ]
[0.276 0.07 1.108 0.683 0. 0. ]
[0.133 0.684 0.113 0.14 0. 0. ]
[0.996 0.211 1.97 2.003 0.449 1.055]
[1.297 1.69 0.73 1.478 0.497 0.586]]
Iteration 200 Node 5
[[1.296 1.345 2.085 1.7 1.057 0.355]
[0.277 0.069 1.109 0.68 0. 0. ]
[0.133 0.69 0.105 0.132 0. 0. ]
[0.996 0.203 1.977 2.01 0.448 1.06 ]
[1.297 1.694 0.724 1.479 0.495 0.585]]
v [array([[ 0.391, 0.46 , 1.711, 1.379, 0.42 , -0.219],
[-0.189, -0.372, 0.936, 0.547, -0.305, -0.284],
[ 0.591, 1.109, 0.311, 0.319, 0.31 , 0.292],
[ 1.454, 0.683, 2.119, 2.13 , 0.761, 1.329],
[ 1.754, 2.121, 0.923, 1.626, 0.814, 0.882]]), array([[ 1.256, 1.247, 1.37 , 1.349, 0.707, 0.686],
[ 0.343, 0.358, 1.028, 1.046, 0.089, 0.113],
[-0.113, -0.096, 0.873, 0.904, -0.2 , -0.182],
[ 1.257, 1.24 , 1.371, 1.353, 0.703, 0.695],
[ 1.257, 1.251, 1.358, 1.348, 0.702, 0.689]]), array([[ 1.347, 1.332, 1.407, 1.384, 0.769, 0.752],
[ 1.35 , 1.332, 1.402, 1.378, 0.766, 0.756],
[-0.022, -0.01 , 0.909, 0.939, -0.137, -0.115],
[-0.022, 0.01 , 0.888, 0.917, -0.161, -0.148],
[ 1.348, 1.337, 1.394, 1.382, 0.764, 0.755]]), array([[ 1.255, 1.241, 1.375, 1.358, 0.711, 0.697],
[ 1.258, 1.241, 1.37 , 1.352, 0.709, 0.701],
[ 1.257, 1.248, 1.361, 1.347, 0.708, 0.701],
[-0.113, -0.081, 0.856, 0.891, -0.219, -0.203],
[ 0.343, 0.35 , 1.038, 1.052, 0.091, 0.104]]), array([[-0.112, -0.087, 0.857, 0.895, -0.222, -0.193],
[ 1.258, 1.245, 1.37 , 1.35 , 0.71 , 0.696],
[ 1.256, 1.251, 1.361, 1.345, 0.71 , 0.696],
[ 1.256, 1.238, 1.375, 1.36 , 0.71 , 0.701],
[ 0.342, 0.353, 1.037, 1.05 , 0.092, 0.1 ]]), array([[ 0.663, 0.608, 0.481, 0.835, 0.015, 0.678],
[ 0.781, 0.995, 1.094, 1.527, 0.432, 0.419],
[ 1.832, 1.298, 2.384, 2.346, 1.01 , 1.008],
[ 0.969, 1.71 , 0.591, 0.549, 0.607, 0.025],
[-0.245, -0.612, 1.45 , 0.743, -0.063, -0.13 ]])]
Iteration 300 Node 0
[[1.267 1.288 2.301 1.77 1.186 0.187]
[0.377 0.017 1.122 0.503 0. 0. ]
[0.071 0.845 0. 0. 0. 0. ]
[0.95 0. 2.189 2.249 0.39 1.263]
[1.335 1.85 0.387 1.477 0.424 0.55 ]]
Iteration 300 Node 1
[[1.267 1.286 2.306 1.771 1.189 0.183]
[0.381 0.015 1.122 0.499 0. 0. ]
[0.07 0.847 0. 0. 0. 0. ]
[0.947 0. 2.193 2.253 0.388 1.266]
[1.336 1.853 0.379 1.477 0.423 0.55 ]]
Iteration 300 Node 2
[[1.267 1.283 2.31 1.772 1.193 0.18 ]
[0.384 0.013 1.123 0.495 0. 0. ]
[0.068 0.849 0. 0. 0. 0. ]
[0.944 0. 2.196 2.257 0.386 1.27 ]
[1.338 1.855 0.371 1.476 0.422 0.55 ]]
Iteration 300 Node 3
[[1.267 1.281 2.314 1.773 1.196 0.177]
[0.387 0.011 1.123 0.49 0. 0. ]
[0.066 0.85 0. 0. 0. 0. ]
[0.941 0. 2.2 2.261 0.383 1.273]
[1.339 1.857 0.363 1.476 0.421 0.55 ]]
Iteration 300 Node 4
[[1.267 1.279 2.319 1.774 1.2 0.174]
[0.39 0.009 1.123 0.486 0. 0. ]
[0.065 0.852 0. 0. 0. 0. ]
[0.938 0. 2.203 2.265 0.381 1.276]
[1.341 1.86 0.355 1.475 0.42 0.55 ]]
Iteration 300 Node 5
[[1.267 1.276 2.323 1.775 1.203 0.171]
[0.393 0.008 1.123 0.482 0. 0. ]
[0.063 0.854 0. 0. 0. 0. ]
[0.935 0. 2.206 2.269 0.378 1.28 ]
[1.342 1.862 0.348 1.475 0.419 0.549]]
v [array([[ 0.357, 0.399, 1.957, 1.462, 0.563, -0.408],
[-0.081, -0.429, 0.956, 0.355, -0.305, -0.284],
[ 0.527, 1.29 , 0.171, 0.152, 0.309, 0.292],
[ 1.405, 0.444, 2.362, 2.403, 0.699, 1.557],
[ 1.792, 2.296, 0.554, 1.629, 0.733, 0.843]]), array([[ 1.256, 1.244, 1.368, 1.347, 0.707, 0.686],
[ 0.342, 0.356, 1.026, 1.044, 0.089, 0.113],
[-0.109, -0.094, 0.882, 0.913, -0.2 , -0.182],
[ 1.253, 1.245, 1.367, 1.35 , 0.701, 0.693],
[ 1.257, 1.249, 1.356, 1.346, 0.703, 0.689]]), array([[ 1.345, 1.329, 1.404, 1.381, 0.768, 0.751],
[ 1.348, 1.329, 1.399, 1.376, 0.764, 0.754],
[-0.02 , -0.009, 0.918, 0.947, -0.136, -0.114],
[-0.018, 0.018, 0.888, 0.916, -0.158, -0.144],
[ 1.346, 1.333, 1.391, 1.38 , 0.763, 0.754]]), array([[ 1.255, 1.24 , 1.372, 1.355, 0.711, 0.696],
[ 1.258, 1.24 , 1.368, 1.35 , 0.707, 0.699],
[ 1.253, 1.244, 1.368, 1.354, 0.707, 0.699],
[-0.108, -0.072, 0.856, 0.891, -0.215, -0.199],
[ 0.342, 0.348, 1.035, 1.05 , 0.091, 0.104]]), array([[-0.111, -0.088, 0.855, 0.894, -0.221, -0.192],
[ 1.259, 1.244, 1.369, 1.349, 0.71 , 0.696],
[ 1.255, 1.247, 1.368, 1.353, 0.71 , 0.696],
[ 1.253, 1.245, 1.372, 1.357, 0.708, 0.7 ],
[ 0.344, 0.352, 1.036, 1.048, 0.093, 0.101]]), array([[ 0.699, 0.677, 0.243, 0.761, -0.127, 0.867],
[ 0.674, 1.06 , 1.082, 1.727, 0.435, 0.422],
[ 1.894, 1.122, 2.494, 2.481, 1.01 , 1.009],
[ 1.015, 1.92 , 0.354, 0.283, 0.665, -0.206],
[-0.282, -0.779, 1.827, 0.748, 0.017, -0.091]])]
Iteration 400 Node 0
[[1.28 1.148 2.529 1.828 1.362 0.031]
[0.485 0. 1.106 0.271 0. 0. ]
[0.013 0.906 0. 0. 0. 0. ]
[0.801 0. 2.359 2.446 0.262 1.425]
[1.422 1.946 0.005 1.456 0.376 0.544]]
Iteration 400 Node 1
[[1.281 1.145 2.532 1.829 1.365 0.028]
[0.487 0. 1.104 0.266 0. 0. ]
[0.012 0.907 0. 0. 0. 0. ]
[0.798 0. 2.361 2.45 0.259 1.429]
[1.423 1.948 0.002 1.455 0.375 0.543]]
Iteration 400 Node 2
[[1.281 1.142 2.535 1.83 1.369 0.025]
[0.489 0. 1.102 0.261 0. 0. ]
[0.011 0.909 0. 0. 0. 0. ]
[0.795 0. 2.362 2.454 0.257 1.432]
[1.424 1.949 0.001 1.454 0.374 0.543]]
Iteration 400 Node 3
[[1.282 1.14 2.537 1.832 1.373 0.022]
[0.491 0. 1.099 0.257 0. 0. ]
[0.01 0.91 0. 0. 0. 0. ]
[0.792 0. 2.363 2.458 0.255 1.435]
[1.425 1.95 0. 1.454 0.372 0.542]]
Iteration 400 Node 4
[[1.282 1.137 2.539 1.833 1.377 0.019]
[0.493 0. 1.096 0.252 0. 0. ]
[0.009 0.911 0. 0. 0. 0. ]
[0.79 0. 2.364 2.462 0.252 1.439]
[1.427 1.951 0. 1.453 0.371 0.542]]
Iteration 400 Node 5
[[1.282 1.134 2.541 1.834 1.38 0.016]
[0.495 0. 1.093 0.247 0. 0. ]
[0.008 0.913 0. 0. 0. 0. ]
[0.787 0. 2.365 2.466 0.25 1.442]
[1.428 1.953 0. 1.453 0.37 0.542]]
v [array([[ 0.37 , 0.259, 2.184, 1.52 , 0.739, -0.563],
[ 0.029, -0.444, 0.94 , 0.123, -0.304, -0.284],
[ 0.468, 1.351, 0.171, 0.151, 0.309, 0.292],
[ 1.255, 0.444, 2.531, 2.599, 0.57 , 1.719],
[ 1.878, 2.391, 0.175, 1.607, 0.686, 0.836]]), array([[ 1.256, 1.243, 1.366, 1.347, 0.707, 0.686],
[ 0.344, 0.358, 1.025, 1.044, 0.089, 0.113],
[-0.108, -0.094, 0.882, 0.913, -0.2 , -0.181],
[ 1.252, 1.245, 1.365, 1.35 , 0.701, 0.693],
[ 1.256, 1.247, 1.361, 1.345, 0.702, 0.689]]), array([[ 1.345, 1.328, 1.402, 1.381, 0.768, 0.751],
[ 1.346, 1.33 , 1.397, 1.375, 0.763, 0.753],
[-0.019, -0.009, 0.918, 0.947, -0.136, -0.114],
[-0.017, 0.018, 0.886, 0.917, -0.157, -0.143],
[ 1.345, 1.332, 1.398, 1.379, 0.763, 0.753]]), array([[ 1.255, 1.239, 1.37 , 1.356, 0.711, 0.696],
[ 1.256, 1.241, 1.365, 1.35 , 0.707, 0.699],
[ 1.253, 1.243, 1.367, 1.354, 0.706, 0.699],
[-0.107, -0.072, 0.854, 0.891, -0.214, -0.198],
[ 0.343, 0.348, 1.043, 1.05 , 0.091, 0.104]]), array([[-0.11 , -0.089, 0.853, 0.894, -0.221, -0.192],
[ 1.258, 1.245, 1.366, 1.348, 0.71 , 0.696],
[ 1.255, 1.247, 1.368, 1.352, 0.709, 0.695],
[ 1.253, 1.244, 1.369, 1.357, 0.707, 0.699],
[ 0.345, 0.353, 1.044, 1.049, 0.094, 0.101]]), array([[ 0.685, 0.819, 0.025, 0.702, -0.304, 1.021],
[ 0.567, 1.069, 1.108, 1.96 , 0.436, 0.423],
[ 1.952, 1.063, 2.494, 2.482, 1.011, 1.009],
[ 1.163, 1.92 , 0.194, 0.086, 0.793, -0.37 ],
[-0.367, -0.871, 2.179, 0.77 , 0.065, -0.084]])]
Iteration 500 Node 0
[[1.258 1.026 2.618 1.878 1.527 0. ]
[0.604 0. 0.968 0.046 0. 0. ]
[0. 0.931 0. 0. 0. 0. ]
[0.682 0. 2.414 2.657 0.167 1.553]
[1.455 2.043 0. 1.42 0.306 0.447]]
Iteration 500 Node 1
[[1.258 1.024 2.619 1.879 1.53 0. ]
[0.607 0. 0.966 0.041 0. 0. ]
[0. 0.931 0. 0. 0. 0. ]
[0.68 0. 2.415 2.661 0.165 1.555]
[1.456 2.045 0. 1.419 0.305 0.445]]
Iteration 500 Node 2
[[1.257 1.022 2.621 1.88 1.533 0. ]
[0.609 0. 0.963 0.037 0. 0. ]
[0. 0.931 0. 0. 0. 0. ]
[0.678 0. 2.416 2.665 0.163 1.557]
[1.456 2.047 0. 1.418 0.303 0.443]]
Iteration 500 Node 3
[[1.256 1.02 2.623 1.881 1.536 0. ]
[0.612 0. 0.96 0.032 0. 0. ]
[0. 0.931 0. 0. 0. 0. ]
[0.675 0. 2.417 2.67 0.162 1.56 ]
[1.457 2.049 0. 1.417 0.302 0.44 ]]
Iteration 500 Node 4
[[1.255 1.017 2.624 1.881 1.539 0. ]
[0.614 0. 0.957 0.028 0. 0. ]
[0. 0.931 0. 0. 0. 0. ]
[0.673 0. 2.418 2.674 0.16 1.562]
[1.457 2.052 0. 1.417 0.3 0.438]]
Iteration 500 Node 5
[[1.254 1.015 2.626 1.882 1.543 0. ]
[0.617 0. 0.955 0.023 0. 0. ]
[0. 0.931 0. 0. 0. 0. ]
[0.671 0. 2.419 2.678 0.158 1.564]
[1.458 2.054 0. 1.416 0.299 0.436]]
v [array([[ 0.35 , 0.139, 2.273, 1.57 , 0.905, -0.592],
[ 0.147, -0.444, 0.801, -0.103, -0.304, -0.284],
[ 0.456, 1.374, 0.171, 0.151, 0.309, 0.292],
[ 1.137, 0.443, 2.585, 2.81 , 0.475, 1.846],
[ 1.911, 2.487, 0.171, 1.571, 0.615, 0.738]]), array([[ 1.254, 1.242, 1.366, 1.347, 0.707, 0.689],
[ 0.343, 0.357, 1.024, 1.044, 0.089, 0.113],
[-0.106, -0.091, 0.882, 0.913, -0.199, -0.181],
[ 1.253, 1.244, 1.365, 1.35 , 0.702, 0.692],
[ 1.256, 1.247, 1.363, 1.345, 0.702, 0.687]]), array([[ 1.344, 1.327, 1.401, 1.381, 0.767, 0.754],
[ 1.347, 1.329, 1.397, 1.376, 0.763, 0.753],
[-0.018, -0.006, 0.918, 0.948, -0.136, -0.114],
[-0.018, 0.02 , 0.885, 0.916, -0.158, -0.145],
[ 1.345, 1.331, 1.399, 1.379, 0.763, 0.752]]), array([[ 1.254, 1.239, 1.37 , 1.355, 0.71 , 0.699],
[ 1.257, 1.241, 1.365, 1.35 , 0.707, 0.699],
[ 1.254, 1.241, 1.367, 1.354, 0.706, 0.698],
[-0.108, -0.071, 0.854, 0.891, -0.215, -0.2 ],
[ 0.343, 0.35 , 1.044, 1.05 , 0.092, 0.104]]), array([[-0.11 , -0.087, 0.853, 0.894, -0.22 , -0.189],
[ 1.258, 1.245, 1.366, 1.348, 0.709, 0.695],
[ 1.255, 1.245, 1.368, 1.352, 0.709, 0.695],
[ 1.253, 1.244, 1.37 , 1.357, 0.708, 0.698],
[ 0.344, 0.353, 1.044, 1.048, 0.094, 0.1 ]]), array([[ 0.708, 0.94 , -0.063, 0.652, -0.47 , 1.039],
[ 0.449, 1.072, 1.248, 2.184, 0.436, 0.423],
[ 1.96 , 1.036, 2.495, 2.482, 1.011, 1.01 ],
[ 1.283, 1.92 , 0.141, -0.125, 0.888, -0.491],
[-0.4 , -0.968, 2.179, 0.806, 0.135, 0.019]])]
In [ ]:
# Final Algorithm value
V@wta.get_final_surv_prob(q, all_w[0].value)
Out[ ]:
33.779418599372576
In [ ]:
# Final Algorithm solution
all_w[0].value
Out[ ]:
array([[1.258, 1.026, 2.618, 1.878, 1.527, 0. ],
[0.604, 0. , 0.968, 0.046, 0. , 0. ],
[0. , 0.931, 0. , 0. , 0. , 0. ],
[0.682, 0. , 2.414, 2.657, 0.167, 1.553],
[1.455, 2.043, 0. , 1.42 , 0.306, 0.447]])In [ ]:
# Plot log values using log/log scale
import matplotlib.pyplot as plt
plt.plot(log)
plt.yscale("log")
plt.xscale("log")
# Add horizontal line at y = x_full
plt.axhline(y=p, color='red', linestyle='dashed')
plt.xlabel("Iteration")
plt.ylabel("Total Value")
plt.title("Total Value vs. Iteration")
plt.show()
In [ ]:
# Plot change in W values using log/log scale
plt.plot(log_e)
plt.yscale("log")
plt.xscale("log")
# Label axes
plt.xlabel("Iteration")
plt.ylabel("Change in w")
plt.title("Change in w vs. Iteration")
plt.show()
In [ ]: