oldNC_Route.py

# -*- coding: utf-8 -*-
import numpy as np
import cvxpy as cvx

## Initialization

oxide = True
decayHastener = True
maxHeat = 9999

enableTBU = True
enableLE = True
enableHE = True
enableMOX = True

# Concatenate with comma(,) not semicolon(,)
forbiddenFuels = np.array([
    'placeholder'
    ])

# RF / RTG
goal = 'RTG'
print(goal)

# Starting resources

# Ingots
Thorium = 500.
Uranium = 600.

# will use amount specified above in ticks specified below
cycle = 138000.

# Cell multiplier: Cell numbers will be a multiple of this number
multiplier = 1

#will be scaled inverse to the sqrt(heat) of the fuel
cellN = 80

reactorN = 0

efficiencyMin = 100
efficiencyMax = 100

## Tables

TBU = np.zeros((27,1))
LE = np.zeros((27,1))
HE = np.zeros((27,1))
MOX = np.zeros((27,1))

TBU[0] = 1
LE[1:6:2] = 1
LE[9:27:2] = 1
HE[2:7:2] = 1
HE[10:27:2] = 1
MOX[7:9] = 1

fuelOrder = np.array([
    'TBU',      'LEU-235',  'HEU-235',  'LEU-233',  'HEU-233',
    'LEN-236',  'HEN-236',  'MOX-239',  'MOX-241',
    'LEP-239',  'HEP-239',  'LEP-241',  'HEP-241',  'LEA-242',  'HEA-242',
    'LECm-243', 'HECm-243', 'LECm-245', 'HECm-245', 'LECm-247', 'HECm-247',
    'LEB-248',  'HEB-248',  'LECf-249', 'HECf-249', 'LECf-251', 'HECf-251'
    ])

isoOrder = np.array([
    'Th-230','Th-232',
    'U-233','U-235','U-238',
    'Np-236','Np-237',
    'Pu-238','Pu-239','Pu-241','Pu-242',
    'Am-241','Am-242','Am-243',
    'Cm-243','Cm-245','Cm-246','Cm-247',
    'Bk-247','Bk-248','Cf-249','Cf-250','Cf-251','Cf-252'
    ])

heat = np.array([
    18,     50,     300,    60,     360,
    36,     216,    57.5,   97.5,
    40,     240,    70,     420,    94,     564,
    112,    672,    68,     408,    54,     324,
    52,     312,    116,    696,    120,    720
    ])

if oxide == True:
    heat[0:7] *= 1.25
    heat[9:27] *= 1.25

## Simulating reactor efficiency
# eff = -heat;
# eff = -log(heat);
eff = np.min(heat) / np.log(heat) + 1
eff = eff - np.min(eff)
eff = eff / np.max(eff) * (efficiencyMax-efficiencyMin) + efficiencyMin

eff = eff/100

## Taking into account time it takes to burn fuels
tick = np.array([
    144000, 72000,  72000,  64000,  64000,
    102000, 102000, 84000,  56000,
    92000,  92000,  60000,  60000,  54000,  54000,
    52000,  52000,  68000,  68000,  78000,  78000,
    86000,  86000,  60000,  60000,  58000,  58000
    ])

CC = cycle / tick * multiplier

RFtNormal = np.array([[
    60,     120,    480,    144,    576,
    90,     360,    155.4,  243.6,
    105,    420,    165,    660,    192,    768,
    210,    840,    162,    648,    138,    552,
    135,    540,    216,    864,    225,    900]])

RFtNormal = RFtNormal * CC * eff

RFtOxide = RFtNormal*1.4
RFtOxide[0,8] = RFtNormal[0,8]
RFtOxide[0,9] = RFtNormal[0,9]

## Hardcoded arrays
separatorMatrix = np.array([
    [1,9,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
    [0,0,0,1,9,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]])

fuelMatrix = np.array([
    [0,-81,16,8,0,8,32,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
    [0,0,0,-9,-32,0,8,0,8,8,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
    [0,0,0,-36,-25,0,16,0,4,0,24,0,0,0,0,0,0,0,0,0,0,0,0,0],
    [0,0,-9,0,-72,0,0,0,4,4,32,0,0,24,0,0,0,0,0,0,0,0,0,0],
    [0,0,-36,0,-45,32,8,0,0,0,16,0,0,8,0,0,0,0,0,0,0,0,0,0],
    [0,0,0,0,0,-9,-68,0,0,0,32,0,8,20,0,0,0,0,0,0,0,0,0,0],
    [0,0,0,0,16,-36,-45,8,8,0,32,0,0,0,0,0,0,0,0,0,0,0,0,0],
    [0,0,0,0,-32,0,0,0,-9,0,12,0,0,8,4,0,0,0,0,0,0,0,0,0],
    [0,0,0,0,-32,0,0,0,0,-1,8,0,0,0,0,0,8,0,0,0,0,0,0,0],
    [0,0,0,0,0,0,0,0,-1,0,-48,0,0,0,4,0,28,0,0,0,0,0,0,0],
    [0,0,0,0,0,0,0,0,-36,0,-45,8,24,0,0,8,24,0,0,0,0,0,0,0],
    [0,0,0,0,0,0,0,0,0,-9,-68,0,4,8,0,0,48,0,0,0,0,0,0,0],
    [0,0,0,0,0,0,0,0,0,-36,-45,8,0,0,0,8,24,24,0,0,0,0,0,0],
    [0,0,0,0,0,0,0,0,0,0,0,0,-9,-72,8,8,40,8,0,0,0,0,0,0],
    [0,0,0,0,0,0,0,0,0,0,0,0,-36,-45,0,16,32,8,8,0,0,0,0,0],
    [0,0,0,0,0,0,0,0,0,0,0,0,0,0,-9,0,-40,0,16,8,8,0,0,0],
    [0,0,0,0,0,0,0,0,0,0,0,0,0,0,-36,0,-21,0,24,8,8,0,0,0],
    [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-9,-72,0,40,8,4,0,0,12],
    [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-36,-45,0,48,4,4,0,8,0],
    [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-72,-9,20,4,0,0,8,32],
    [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-45,-36,0,8,8,0,24,24],
    [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-72,-9,4,0,4,56],
    [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-45,-36,0,8,8,48],
    [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-9,16,8,-32],
    [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-36,32,16,-29],
    [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-5,-12],
    [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-20,3]])

fuelMatrix = np.diag(CC)@fuelMatrix

decayMatrix = np.array([
    [1,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
    [0,1,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
    [0,0,1,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
    [1,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
    [0,0,0,1,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
    [0,0,0,0,0,0,1,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
    [0,0,0,0,1,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0],
    [0,0,0,0,0,0,1,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0],
    [1,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0],
    [0,0,0,0,0,0,0,0,1,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0],
    [0,0,0,0,0,0,0,0,1,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0],
    [0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,-1,0,0,0,0,0,0,0,0],
    [0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,-1,0,0,0,0,0,0,0],
    [0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,-1,0,0,0,0,0,0],
    [0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,-1,0,0,0,0,0],
    [0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0],
    [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,-1,0,0,0],
    [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,-1,0,0],
    [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,-1,0],
    [0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1]])

#isotopeSeparatorRF = .5 #kRF that takes to process one ingot
#decayHastenerRF = .5 #kRF that takes to process small clumps

A = np.concatenate((separatorMatrix, fuelMatrix, decayMatrix)).T

# if goal == "RF" and oxide == True:
#     f1 = RFtOxide
#     f2 = np.zeros((20))
# elif goal == "RF":
#     f1 = RFtNormal
#     f2 = np.zeros((20))
# elif goal == "RTG":
#     f1 = fuelMatrix[:,21]
#     f2 = np.array([0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0])
#     #HEB-248, LECf-249, HECf-249 produces Cf-250 for these amounts
# f1 = f1[..., np.newaxis]
# f2 = f2[..., np.newaxis]

C0 = heat > maxHeat
maximumN = cellN / np.sqrt(heat/heat[1])

FF = np.squeeze(((enableLE * LE + enableHE * HE + enableMOX * MOX + enableTBU * TBU)==0).T) | \
    np.isin(fuelOrder, forbiddenFuels) | C0

FF = FF[..., np.newaxis]

N = 27
M = N+1
C = np.zeros(49)
Result = np.zeros([50, 27])

x = cvx.Variable((27, 1))
Z = cvx.Variable((27, 1), integer = True)
y = cvx.Variable((20, 1))

tol = 10 ** -6

ub = maximumN / multiplier
ub = ub[..., np.newaxis]

if goal == "RF" and oxide == True:
    f1 = RFtOxide
    objective = cvx.Maximize(f1.T @ x)
elif goal == "RF":
    f1 = RFtNormal
    objective = cvx.Maximize(f1.T @ x)
elif goal == "RTG":
    f1 = fuelMatrix[:,21]
    f2 = np.array([0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0])
    objective = cvx.Maximize(f1.T @ x + f2.T @ y)

constraints = [x >= 0, x <= ub]

if decayHastener == True:
    constraints += [y >= 0]
else:
    constraints += [y == 0]

constraints += [separatorMatrix.T @ np.array([[Thorium], [Uranium]])+
                fuelMatrix.T @ x +
                decayMatrix.T @ y >= 0]

constr_int = [x == Z]

while True:

    FuelConstraints = [cvx.multiply(x, FF) == 0]

    prob = cvx.Problem(objective,
                          constraints + FuelConstraints + constr_int)

    prob.solve()
    N = np.sum(x.value!=0)
    print('Number of reactors: ' + str(N))
    if prob.status == 'infeasible':
        break
    if prob.value == 0:
        break
    if M>N:
        Result[0, N-1] = prob.value
        Result[1, N-1] = Thorium
        Result[2, N-1] = Uranium
        Result[3:30, N-1] = np.squeeze(x.value)
        Result[30:50, N-1] = np.squeeze(y.value)
        M=N

    fArray = np.zeros(27)
    for k in np.where(x.value>tol)[0]:
        FF2 = np.copy(FF)
        FF2[k] = True
        FuelConstraints = [cvx.multiply(x, FF2) == 0]
        prob = cvx.Problem(objective, constraints + FuelConstraints)
        prob.solve()
        fArray[k] = prob.value


    k = np.where(fArray == np.amax(fArray))
    FF[k] = True
    print('Eliminating: ')
    for fuel in fuelOrder[fArray == np.amax(fArray)]:
        print(fuel)