查看“水域自动化笔记ST13血魔法”的源代码

== 超链接 ==
[[水域自动化笔记ST13|返回ST13页面]]

[[水域自动化笔记目录|返回索引页面]]

== 调包 ==
<pre>import numpy as np
from scipy.optimize import minimize
import matplotlib.pyplot as plt
from tqdm import tqdm</pre>

== 参数计算函数 ==
<pre>def get_altar_cap(N_capacity=0, N_aug_capa=0, N_capa_lv2=0, N_aug_capa_lv2=0, base = 10000):
    '''N_capacity=0, N_aug_capa=0, N_capa_lv2=0, N_aug_capa_lv2=0, base = 10000, printing = True'''
    capacity = np.floor(base*(1+N_capacity*0.2+N_capa_lv2*0.4)*(1.075**N_aug_capa)*(1.1556228712619**N_aug_capa_lv2))
    #why aug cap 2 is so strange?
    return capacity

def get_altar_gtpct(N_acc = 0, N_acc_lv2 = 0, IFS_surge_num = 0, default = 20):
    '''N_acc = 0, N_acc_lv2 = 0, IFS_surge_num = 0, default = 20'''
    GTpCT = np.max([1,(default-N_acc-2*N_acc_lv2)])/(1+2*IFS_surge_num)
    return GTpCT

def get_altar_dislocIO(N_disloc = 0, N_disloc_lv2 = 0):
    return np.floor(20*(1.2**N_disloc*1.4**N_disloc_lv2))

def get_altar_dislocIOperGT(N_disloc = 0, N_disloc_lv2 = 0, GTperCT = 20):
    return np.floor(20*(1.2**N_disloc*1.4**N_disloc_lv2)) / GTperCT
    
def get_altar_IO(capacity = 10000, N_disloc = 0, N_disloc_lv2 = 0):
    IO = np.min([capacity/10, np.floor(20*(1.2**N_disloc*1.4**N_disloc_lv2))])
    return IO

def get_altar_IOperGT(capacity = 10000, N_disloc = 0, N_disloc_lv2 = 0, GTperCT = 20):
    IO = np.min([capacity/10, np.floor(20*(1.2**N_disloc*1.4**N_disloc_lv2))])
    return np.min([capacity/10, IO/GTperCT])

def get_altar_maxproc(max_rate = 50*5, N_speed = 0, N_speed_lv2 = 0):
    return max_rate * (1 + 0.2*N_speed + 0.4*N_speed_lv2)
    
def lp_prod_per_op(base = 25, N_sac = 0, N_sac_lv2 = 0):
    return base * (1+0.1*N_sac+0.2*N_sac_lv2)</pre>

== 最优化求解 ==

=== 合成用祭坛 ===

==== 足量强化符文 ====
对于altar_level = 5（五级祭坛），gear_lim = 108（强化符文数，若无则设为0）的祭坛，求解最大化的祭坛充能速率，满足不等式约束条件（>=0）与等式约束条件（==0）

在我的笔记本上N_init=10k时，计算时长约为14min，以下为求解足量强化符文数的代码与初始化。要注意，SLSQP方法是初始值敏感的，错误的初始估计会导致严重偏差的结果

<pre>altar_level = 5
gear_lim = 108
N_init = 10000

rune_slot = [0,4,28,56,108]
rune_limit = rune_slot[altar_level-1]

def target_func(x):
    return - get_altar_maxproc(N_speed = x[8], N_speed_lv2 = x[9])


ineq_cons = {'type': 'ineq',
             'fun' : lambda x: np.array([19-x[4]-2*x[5],
                                         rune_limit - np.sum(x),
                                         gear_lim - (x[2]+x[3]+x[5]+x[7]+x[9]),
                                         get_altar_IOperGT(capacity = get_altar_cap(N_capacity=x[0], N_aug_capa=x[1], 
                                                                                    N_capa_lv2=x[2], N_aug_capa_lv2=x[3]), 
                                                           N_disloc = x[4], N_disloc_lv2 = x[5],
                                                           GTperCT = get_altar_gtpct(N_acc = x[6], N_acc_lv2 = x[7], IFS_surge_num = 2) )
                                        - get_altar_maxproc(N_speed = x[8], N_speed_lv2 = x[9])
                                        ])}

eq_cons = {'type' : 'eq',
           'fun' : lambda x: np.array([np.sum(x) - rune_limit,
                                       x[2]+x[3]+x[5]+x[7]+x[9] - gear_lim,
                                       get_altar_dislocIOperGT(N_disloc = x[4], N_disloc_lv2 = x[5], 
                                                               GTperCT = get_altar_gtpct(N_acc = x[6], N_acc_lv2 = x[7], IFS_surge_num = 2))
                                       - get_altar_cap(N_capacity=x[0], N_aug_capa=x[1], N_capa_lv2=x[2], N_aug_capa_lv2=x[3]) / 10
                                       ])}



best_proc_rate = 0
float_solution = []
ava_count = 0
for i in tqdm(range(N_init), desc='SLSQP'):
    N_capacity = 0+1*np.random.rand()
    N_capa_lv2 = 10+5*np.random.rand()
    N_aug_capa = 0+1*np.random.rand()
    N_aug_capa_lv2 = 2+3*np.random.rand()
    N_disloc = 0+1*np.random.rand()
    N_disloc_lv2 = 10+10*np.random.rand()
    N_acc = 0+2*np.random.rand()
    N_acc_lv2 = 8+1*np.random.rand()
    N_speed = 0+1*np.random.rand()
    N_speed_lv2 = 108 - np.sum([N_capacity, N_aug_capa, N_capa_lv2, N_aug_capa_lv2,
                                N_disloc, N_disloc_lv2, N_acc, N_acc_lv2, N_speed])
    
    if i % 3 == 0 and len(float_solution)!=0:
        x0 = float_solution
    else:
        
        x0 = np.array([N_capacity, N_aug_capa, N_capa_lv2, N_aug_capa_lv2,
                       N_disloc, N_disloc_lv2, N_acc, N_acc_lv2, N_speed, N_speed_lv2])
    
    bounds = [(0, 108) for _ in range(len(x0))]  
    
    sol = minimize(target_func, x0, method='SLSQP', 
                   constraints=[eq_cons, ineq_cons], options={'maxiter': 10000, 'ftol': 1e-9, 'disp': False},
                   bounds=bounds)

    x = sol.x
    IO = get_altar_IOperGT(capacity = get_altar_cap(N_capacity=x[0], N_aug_capa=x[1],N_capa_lv2=x[2], N_aug_capa_lv2=x[3]), 
                           N_disloc = x[4], N_disloc_lv2 = x[5],
                           GTperCT = get_altar_gtpct(N_acc = x[6], N_acc_lv2 = x[7], IFS_surge_num = 2) )
    ava_solution = (-sol.fun <= IO)
    
    if -sol.fun > best_proc_rate and ava_solution:
        ava_count+=1
        float_solution = sol.x
        best_proc_rate = -sol.fun
        sol_opt = sol</pre>

==== 0强化符文 ====
对于无强化符文的最优解搜寻，需要将gear_lim设为0.1，使用如下初始化，运行时间减少很多

<pre>    N_capacity = 7+2*np.random.rand()
    N_capa_lv2 = 0
    N_aug_capa = 0.5+1*np.random.rand()
    N_aug_capa_lv2 = 0
    N_disloc = 24+8*np.random.rand()
    N_disloc_lv2 = 0
    N_acc = 12+8*np.random.rand()
    N_acc_lv2 = 0
    N_speed = 108-N_capacity-N_aug_capa-N_disloc-N_acc
    N_speed_lv2 = 0</pre>

=== 连续结果显示 ===

==== 代码格 ====

<pre>x = sol_opt.x
IO = get_altar_IOperGT(capacity = get_altar_cap(N_capacity=x[0], N_aug_capa=x[1],N_capa_lv2=x[2], N_aug_capa_lv2=x[3]), 
                       N_disloc = x[4], N_disloc_lv2 = x[5],
                       GTperCT = get_altar_gtpct(N_acc = x[6], N_acc_lv2 = x[7], IFS_surge_num = 2) )

print(f"available rate: {ava_count / N_init*100} %")
print("Optimal Variables: ", sol_opt.x)
print("io: ", IO)
print("Maximum Value of the Objective Function: ", -sol_opt.fun)
print(sum(x))</pre>

==== 样例输出 ====

<pre>available rate: 0.1 %
Optimal Variables:  [ 0.38715085  0.30000646 10.48624571  2.16824612  0.43213499 12.72563799
  1.91935781  8.82011957  0.92436075 69.83673975]
io:  7372.1
Maximum Value of the Objective Function:  7279.892012605729
108.0</pre>

=== 手动修改 ===

要将SLSQP返回的浮点数结果转换为整数结果，可以写一格自动化代码，但是也可以手改，手动修改数组的元素运行代码看是否满足条件即可

==== 代码格 ====

<pre>x = [0,0, 4, 8, 0, 13, 1, 9, 0, 73]
print(sum(x))
print('cap', get_altar_cap(N_capacity=x[0], N_aug_capa=x[1],N_capa_lv2=x[2], N_aug_capa_lv2=x[3]))
print('io', get_altar_IOperGT(capacity = get_altar_cap(N_capacity=x[0], N_aug_capa=x[1],N_capa_lv2=x[2], N_aug_capa_lv2=x[3]), 
                       N_disloc = x[4], N_disloc_lv2 = x[5],
                       GTperCT = get_altar_gtpct(N_acc = x[6], N_acc_lv2 = x[7], IFS_surge_num = 2) ))

print('proc rate', get_altar_maxproc(N_speed = x[8], N_speed_lv2 = x[9]))</pre>

==== 样例输出 ====

<pre>108
cap 82699.0
io 7935.0
proc rate 7550.000000000001</pre>

=== 源质生产用祭坛 ===

==== 基于目标IO速率，最小化符文用量 ====
0强化符文版本，运行仅需半分钟

<pre>altar_level = 5
target_IO_rate = 2500
gear_lim = 0.1
N_init = 10000

rune_slot = [0,4,28,56,108]
rune_limit = rune_slot[altar_level-1]

def target_func(x):
    return sum(x)


ineq_cons = {'type': 'ineq',
             'fun' : lambda x: np.array([19-x[4]-2*x[5],
                                         rune_limit - np.sum(x),
                                         gear_lim - (x[2]+x[3]+x[5]+x[7]),
                                         
                                        ])}

eq_cons = {'type' : 'eq',
           'fun' : lambda x: np.array([gear_lim - (x[2]+x[3]+x[5]+x[7]),
                                       get_altar_dislocIOperGT(N_disloc = x[4], N_disloc_lv2 = x[5], 
                                                               GTperCT = get_altar_gtpct(N_acc = x[6], N_acc_lv2 = x[7], IFS_surge_num = 2))
                                       - get_altar_cap(N_capacity=x[0], N_aug_capa=x[1], N_capa_lv2=x[2], N_aug_capa_lv2=x[3]) / 10,

                                       get_altar_IOperGT(capacity = get_altar_cap(N_capacity=x[0], N_aug_capa=x[1], 
                                                                                    N_capa_lv2=x[2], N_aug_capa_lv2=x[3]), 
                                                           N_disloc = x[4], N_disloc_lv2 = x[5],
                                                           GTperCT = get_altar_gtpct(N_acc = x[6], N_acc_lv2 = x[7], IFS_surge_num = 2) )
                                        - target_IO_rate
                                       ])}



minimum_rune_usage = np.inf
float_solution = []
ava_count = 0
for i in tqdm(range(N_init), desc='SLSQP'):
    N_capacity = 7+4*np.random.rand()
    N_capa_lv2 = 0
    N_aug_capa = 1+1*np.random.rand()
    N_aug_capa_lv2 = 0
    N_disloc = 16+8*np.random.rand()
    N_disloc_lv2 = 0
    N_acc = 15+6*np.random.rand()
    N_acc_lv2 = 0
    
    if i % 3 == 0 and len(float_solution)!=0:
        x0 = float_solution
    else:
        
        x0 = np.array([N_capacity, N_aug_capa, N_capa_lv2, N_aug_capa_lv2,
                       N_disloc, N_disloc_lv2, N_acc, N_acc_lv2])
    
    bounds = [(0, 108) for _ in range(len(x0))]  
    
    sol = minimize(target_func, x0, method='SLSQP', 
                   constraints=[eq_cons, ineq_cons], options={'maxiter': 10000, 'ftol': 1e-9, 'disp': False},
                   bounds=bounds)

    x = sol.x
    IO = get_altar_IOperGT(capacity = get_altar_cap(N_capacity=x[0], N_aug_capa=x[1],N_capa_lv2=x[2], N_aug_capa_lv2=x[3]), 
                           N_disloc = x[4], N_disloc_lv2 = x[5],
                           GTperCT = get_altar_gtpct(N_acc = x[6], N_acc_lv2 = x[7], IFS_surge_num = 2) )
    
    ava_solution = (IO >= target_IO_rate) and (sum(x) <= rune_limit+0.5)
    
    if sol.fun < minimum_rune_usage and ava_solution:
        ava_count+=1
        float_solution = sol.x
        best_proc_rate = -sol.fun
        sol_opt = sol
</pre>
=== 连续结果显示 ===

==== 代码格 ====
<pre>print(f"available rate: {ava_count / N_init*100} %")

x = sol_opt.x
IO = get_altar_IOperGT(capacity = get_altar_cap(N_capacity=x[0], N_aug_capa=x[1],N_capa_lv2=x[2], N_aug_capa_lv2=x[3]), 
                       N_disloc = x[4], N_disloc_lv2 = x[5],
                       GTperCT = get_altar_gtpct(N_acc = x[6], N_acc_lv2 = x[7], IFS_surge_num = 2) )


print("Optimal Variables: ", sol_opt.x)
print("io: ", IO)
print("Maximum Value of the Objective Function: ", sol_opt.fun)
print(sum(x))</pre>

==== 样例输出 ====

<pre>available rate: 62.53999999999999 %
Optimal Variables:  [ 7.73967589  1.33906908  0.          0.         23.82983743  0.
 17.35302259  0.        ]
io:  2807.0
Maximum Value of the Objective Function:  50.26160499370298
50.26160499370298</pre>

=== 手动修改 ===

修改x参数可以看到剩余符文、强化符文数和并行数，一般建议并行数控制在150以下

==== 代码格 ====

<pre>x = [8, 0, 0,0, 18, 0, 19, 0]


rune_slot = [0,4,28,56,108]
rune_limit = rune_slot[altar_level-1]

cap = get_altar_cap(N_capacity=x[0], N_aug_capa=x[1],N_capa_lv2=x[2], N_aug_capa_lv2=x[3])
IO = get_altar_IOperGT(capacity = cap, 
                       N_disloc = x[4], N_disloc_lv2 = x[5],
                       GTperCT = get_altar_gtpct(N_acc = x[6], N_acc_lv2 = x[7], IFS_surge_num = 2) )

print(f"cap = {cap} mb, io = {IO} mb/t")
print(f"sum of rune = {sum(x)}, remaining {rune_limit - sum(x)}")
print(f"remaining gear = {np.floor(gear_lim - (x[2]+x[3]+x[5]+x[7]))}")
N_sac_lv2 = np.min([rune_limit - sum(x), np.floor(gear_lim - (x[2]+x[3]+x[5]+x[7]))])
N_sac = rune_limit - sum(x) - N_sac_lv2
print(f"with sac = {N_sac} and sac lv2 = {N_sac_lv2}, lp per op per entity = {lp_prod_per_op(N_sac = N_sac, N_sac_lv2 = N_sac_lv2)} mb")
print(f"the required parallel num for well of suf = {np.ceil(IO*10 / lp_prod_per_op(N_sac = N_sac, N_sac_lv2 = N_sac_lv2))}")</pre>

==== 样例输出 ====

<pre>cap = 26000.0 mb, io = 2600.0 mb/t
sum of rune = 45, remaining 63
remaining gear = 0.0
with sac = 63.0 and sac lv2 = 0.0, lp per op per entity = 182.50000000000003 mb
the required parallel num for well of suf = 143.0</pre>