这是完整代码import math import random import numpy as np import matplotlib.pyplot as plt #import self as self epsilon = 0.5 gamma = 0.1 lr = 0.1 zeros_vector=[] x = []; y = []; X = []; Y = []; agent=[x,y]; object=[X,Y]; random.seed(70) for i in range(10): x.append(random.uniform(0, 1)) y.append(random.uniform(0, 1)) X.append(random.uniform(1, 10)) Y.append(random.uniform(1, 10)) distance = [] for i in range(len(agent[0])): distance_vector = [] for j in range(len(object[0])): dx = agent[0][i] - object[0][j] dy = agent[1][i] - object[1][j] distance_vector.append(math.sqrt(dx * dx + dy * dy)) distance.append(distance_vector) R_table = np.zeros((10, 10)) for i in range(len(agent[0])): for j in range(len(object[0])): R_table[i,j] = 20-distance[i][j] space = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9] #Q_table = [] Q_table = np.zeros((10, 10)) # 进行训练同时测试训练成果 iterate_results = [] # 保存每次测试结果 for i in range(500): print(f"开始第{i + 1}回合。。。") # 初始位置 path = [] # 每个回合要获取10个位置 for j in range(10): remain = set(space) - set(path) # 剩余节点 # s = path[0] # 当前位置 # s_row = Q_table[s] # 当前位置对应的Q表中的行 max_value = -1000 # 在剩余动作中遍历最大值 for rm in remain: Q = Q_table[j][rm] if Q > max_value: max_value = Q a = rm # 随机选择下一行动 if np.random.uniform() < epsilon: a = np.random.choice(np.array(list(set(space) - set(path)))) # 更新Q表 if j != 10: Q_table[j][a] = (1 - lr) * Q_table[j][a] + lr * (R_table[j][a] + gamma * max_value) else: Q_table[j][a] = (1 - lr) * Q_table[j][a] + lr * R_table[j][a] path.append(a) # print(Q_table[j][a]) # 根据当前Q表获取最佳路径距离 result = [] lengths=[] final_length=1000 final_result=[] for k in range(10): loc = k remain1 = set(space) - set(result) # 剩余节点 max_value1 = -1000 # 在剩余节点中遍历最大值 # a = 0 # 利用贪婪策略选择下一行动 for v in remain1: Q_ = Q_table[loc][v] if Q_ > max_value1: a = v max_value1 = Q_ result.append(a) length = 0 for v in range(1, 10): length += distance[result[v - 1]][v-1] print(f"距离为：{length}") if length < final_length: final_length=length final_result=result iterate_results.append(length)

根据您提供的代码，if条件句应该是可以生效的。代码中有一个变量叫做`final_length`，它的初始值为1000。在代码的最后，有一个循环用于计算路径长度，并与`final_length`进行比较，如果路径长度小于`final_length`，则将`final_length`更新为当前路径长度，并将路径保存为`final_result`。所以，如果路径长度小于`final_length`，if条件句应该会被执行。

具体来说，下面这段代码就是if条件句：

if length < final_length:
    final_length=length
    final_result=result

这段代码的作用是，如果当前路径长度`length`小于`final_length`，则将当前路径长度赋值给`final_length`，并将当前路径保存为`final_result`。

如果您有任何其他问题，请随时问我。

相关问题

import numpy as np import matplotlib.pyplot as plt import sympy from scipy.interpolate import interp1d gamma = 1.2 R = 8.314 T0 = 500 Q = 50 * R * T0 a0 = np.sqrt(gamma * R * T0) M0 = 6.216 P_P0 = sympy.symbols('P_P0') num = 81 x0 = np.linspace(0,1,num) t_t0 = np.linspace(0,15,num) x = x0[1:] T_T0 = t_t0[1:] h0 = [] h1 = []#创建拉姆达为1的空数组 r = [] t = [] c = [] s = [] i = 0 for V_V0 in x: n1 = sympy.solve(1 / (gamma-1) * (P_P0 * V_V0 - 1) - 0.5 * (P_P0 + 1) * (1 - V_V0)- gamma * 0 * Q / a0 ** 2,P_P0)#lamuda=0的Hugoniot曲线方程 n2 = sympy.solve(1 / (gamma-1) * (P_P0 * V_V0 - 1) - 0.5 * (P_P0 + 1) * (1 - V_V0)- gamma * 1 * Q / a0 ** 2,P_P0)#lamuda=1的Hugoniot曲线方程 n3 = sympy.solve(-1 * P_P0 + 1 - gamma * M0 ** 2 * (V_V0 - 1),P_P0)#Reyleigh曲线方程 n4 = 12.014556 / V_V0#等温线 n5 = sympy.solve((P_P0 - 1 / (gamma+1) )* (V_V0-gamma / (gamma + 1)) - gamma / ((gamma + 1) ** 2),P_P0)#声速线 n6 = 10.6677 / np.power(V_V0,1.2)#等熵线 h0.append(n1) h1.append(n2) r.append(n3) t.append(n4) c.append(n5) s.append(n6) i = i+1 h0 = np.array(h0) h1 = np.array(h1) r = np.array(r) t = np.array(t) c = np.array(c) s = np.array(s) plt.plot(x,r,label='Rayleigh') plt.plot(x,t,color='purple',label='isothermal') plt.plot(x,s,color='skyblue',label='isentropic') a = np.where(h0 < 0) b = np.where(c < 0) h0 = np.delete(h0,np.where(h0 < 0)[0],axis = 0)#去除解小于0的值 h1 = np.delete(h1,np.where(h1 < 0)[0],axis = 0)#去除解小于0的值 c = np.delete(c,np.where(c < 0)[0],axis = 0)#去除解小于0的值 x0 = np.delete(x,a,axis = 0)#对应去除x轴上错误值的坐标 x1 = np.delete(x,b,axis = 0) plt.plot(x0,h0,label='Hugoniot(lambda=0)') plt.plot(x0,h1,label='Hugoniot(lambda=1)') plt.plot(x1,c,color='yellow',label='soniclocus') plt.ylim((0,50)) plt.legend() # 显示图例 plt.xlabel('V/V0') plt.ylabel('P/P0') f1 = interp1d(x1, c.T, kind='cubic') f2 = interp1d(x,r.T,kind='cubic') f3 = interp1d(x, t.T, kind='cubic') epsilon = 0.0001 x0 = 0.56 y0 = f1(x0) - f2(x0) while abs(y0) > epsilon: df = (f1(x0 + epsilon) - f2(x0 + epsilon) - y0) / epsilon x0 -= y0 / df y0 = f1(x0) - f2(x0) plt.scatter(x0, y0, 50, color ='red') plt.show()

这是一个 Python 代码，我可以解读出来。这段代码主要是用来绘制气体动力学中的一些图形，包括 Hugoniot 曲线、Rayleigh 曲线、等温线、等熵线和声速线等等。其中用到了一些科学计算库，比如 numpy、matplotlib、sympy 和 scipy.interpolate。这段代码还包括了一些数据处理的操作，比如删除数组中小于0的值、插值等等。最后，这段代码还用牛顿迭代法求解了两条曲线的交点。

这个完整的MATLAB代码，可以用python来实现吗，怎么把这个代码转换成python代码，具体代码为？请详细写出

将一个复杂的MATLAB代码转换为Python代码需要逐行翻译和调整语法。以下是该MATLAB代码的主要部分转换为Python的示例：

### Python 代码转换

import numpy as np
from scipy.stats import uniform
import matplotlib.pyplot as plt

# 清理工作空间
np.random.seed(0)

# 问题定义
def myfun(x):
    # 这里假设 `myfun` 已经定义好
    return np.sum(x**2), x

WPc = 16  # 建筑废料生成地点数量
RCc = 6   # 回收中心数量
MMc = 4   # 再生厂商数量
nVar = RCc + RCc * WPc + MMc + MMc * RCc  # 决策变量个数
VarSize = [1, nVar]  # 决策变量矩阵大小
VarMin = 0  # 变量下界
VarMax = 1  # 变量上界
VelMax = (VarMax - VarMin) / 5  # 粒子速度步长极限

# MOPSO 设置
MaxIt = 200  # 最大迭代次数
nPop = 100   # 种群大小
nRep = 100   # 存储规模
w = 0.5      # 惯性权重
wdamp = 0.99 # 惯性权重衰减率
c1 = 1       # 个人学习系数
c2 = 2       # 全局学习系数
nGrid = 7    # 每个维度的栅格数
alpha = 0.1  # 膨胀系数
beta = 2     # 领导者选择压力
gamma = 2    # 删除选择压力
mu = 0.1     # 变异率

# 初始化
wmax = 0.9
wmin = 0.4

# 混沌部分
u0 = 0.5
y0 = 0.5
U_s = [u0]
for _ in range(2 * MaxIt - 1):
    y = np.cos(np.pi * 2 * u0) + y0 * np.exp(-3)
    u = (u0 + 400 + 12 * y) - np.floor(u0 + 400 + 12 * y)
    u0 = u
    y0 = y
    U_s.append(u0)

class Particle:
    def __init__(self):
        self.Velocity = uniform.rvs(loc=VarMin, scale=VarMax - VarMin, size=VarSize)
        self.Position = uniform.rvs(loc=VarMin, scale=VarMax - VarMin, size=VarSize)
        self.Best = {'Position': self.Position.copy(), 'Cost': None}
        self.Cost = None

particles = [Particle() for _ in range(nPop)]

for i in range(nPop):
    fit, pp = myfun(particles[i].Position)
    particles[i].Position = pp
    particles[i].Cost = fit
    particles[i].Best['Position'] = particles[i].Position.copy()
    particles[i].Best['Cost'] = particles[i].Cost

def DetermineDomination(particles):
    for p in particles:
        p.Dominated = False
        for q in particles:
            if all(p.Cost <= q.Cost) and any(p.Cost < q.Cost):
                p.Dominated = True
                break
    return [p for p in particles if not p.Dominated]

particles = DetermineDomination(particles)
rep = GetNonDominatedParticles(particles)
rep_costs = GetCosts(rep)
G = CreateHypercubes(rep_costs, nGrid, alpha)

for r in rep:
    r.GridIndex, r.GridSubIndex = GetGridIndex(r, G)

# MOPSO 主循环
for it in range(1, MaxIt + 1):
    w = wmax - ((wmax - wmin) / MaxIt) * it
    for i in range(nPop):
        rep_h = SelectLeader(rep, beta)
        index = np.random.permutation(2 * MaxIt)[:2]
        particles[i].Velocity = w * particles[i].Velocity \
                                + c1 * U_s[index[0]] * (particles[i].Best['Position'] - particles[i].Position) \
                                + c2 * U_s[index[1]] * (rep_h.Position - particles[i].Position)
        particles[i].Velocity = np.clip(particles[i].Velocity, -VelMax, VelMax)
        particles[i].Position += particles[i].Velocity
        flag = (particles[i].Position < VarMin) | (particles[i].Position > VarMax)
        particles[i].Velocity[flag] = -particles[i].Velocity[flag]
        particles[i].Position = np.clip(particles[i].Position, VarMin, VarMax)
        
        fit, pp = myfun(particles[i].Position)
        particles[i].Position = pp
        particles[i].Cost = fit
        
        # 变异操作
        if np.random.rand() < mu:
            flag = np.random.randint(1, 4)
            NewSol = operator_3(particles[i].Position, 1, nVar, flag)
            fit, pp = myfun(NewSol)
            NewSol = pp
            NewSol_Cost = fit
            
            if Dominates({'Position': NewSol, 'Cost': NewSol_Cost}, particles[i]):
                particles[i].Position = NewSol
                particles[i].Cost = NewSol_Cost
            elif Dominates(particles[i], {'Position': NewSol, 'Cost': NewSol_Cost}):
                pass
            else:
                if np.random.rand() < 0.5:
                    particles[i].Position = NewSol
                    particles[i].Cost = NewSol_Cost
        
        if Dominates(particles[i], particles[i].Best):
            particles[i].Best['Position'] = particles[i].Position.copy()
            particles[i].Best['Cost'] = particles[i].Cost
        elif not Dominates(particles[i].Best, particles[i]):
            if np.random.rand() < 0.5:
                particles[i].Best['Position'] = particles[i].Position.copy()
                particles[i].Best['Cost'] = particles[i].Cost
    
    particles = DetermineDomination(particles)
    nd_particle = GetNonDominatedParticles(particles)
    rep.extend(nd_particle)
    rep = DetermineDomination(rep)
    rep = GetNonDominatedParticles(rep)
    
    for r in rep:
        r.GridIndex, r.GridSubIndex = GetGridIndex(r, G)
    
    if len(rep) > nRep:
        EXTRA = len(rep) - nRep
        rep = DeleteFromRep(rep, EXTRA, gamma)
        rep_costs = GetCosts(rep)
        G = CreateHypercubes(rep_costs, nGrid, alpha)
    
    print(f'IMOPSO Iteration {it}: Number of Repository Particles = {len(rep)}')

# 结果
costs = GetCosts(particles)
rep_costs = GetCosts(rep)
plt.plot(rep_costs[:, 1], rep_costs[:, 0], 'ro')
plt.grid(True)
plt.show()

# 实验数据
WP = np.array([
    [22, 82, 1043.061012], [49, 86, 2005.26846], [17, 55, 5897.871275],
    [22, 77, 1770.120257], [24, 55, 9544.653735], [45, 67, 9074.48371],
    [55, 82, 2220.08252], [69, 81, 7337.41651], [26, 36, 1377.185126],
    [38, 51, 3514.933247], [37, 44, 2154.885524], [59, 44, 3631.591192],
    [62, 44, 1603.743041], [74, 50, 12699.22979], [33, 24, 310.9428793],
    [38, 23, 1522.232645]
])

RC = np.array([
    [62, 95, 6000, 20000000, 65, 7.4, 0.81], [62, 34, 8000, 35000000, 80, 8.5, 0.88],
    [16, 30, 7700, 22000000, 70, 7.8, 0.78], [60, 66, 6500, 26000000, 75, 7.1, 0.821],
    [68, 7500, 22000000, 70, 7.3, 0.80], [71, 22, 7000, 22000000, 75, 7.4, 0.82]
])

MM = np.array([
    [46, 99, 2800, 29000000, 180, 56, 0.80], [57, 55, 3100, 31000000, 185, 52, 0.82],
    [58, 24, 3500, 33000000, 190, 50, 0.85], [6, 45, 4000, 35000000, 195, 60, 0.85]
])

SM = np.array([20, 52, 75, 1.53])  # 二级市场参数
LF = np.array([75, 73, 20, 3.58])  # 填埋区的参数
S = 6000  # 公路运输能力
alpha = 3.5  # 单位废料运输成本
beta = 0.228  # 单位建筑废料单位距离运输碳排放因子
q1 = 1
q2 = 4
q3 = 1
q4 = 2
epsilon_SM = 0.2
epsilon_MM = 0.4
T = 24

X = np.random.rand(1, RCc + RCc * WPc + MMc + MMc * RCc)

# 选择回收中心
chosen_RC = X[0, :RCc]
code = chosen_RC.copy()
chosen_RC[chosen_RC > 0.5] = 1
chosen_RC[chosen_RC == 0.5] = 1
chosen_RC[chosen_RC < 0.5] = 0

while np.sum(chosen_RC) < q1:
    index = np.where(chosen_RC == 0)[0]
    guodu = len(index)
    guoduin = np.random.randint(0, guodu)
    chosen_RC[index[guoduin]] = 1
    guodu = code[index[guoduin]] + 0.5
    if guodu > 1:
        code[index[guoduin]] = 1
    else:
        code[index[guoduin]] = guodu

while np.sum(chosen_RC) > q2:
    index = np.where(chosen_RC == 1)[0]
    guodu = len(index)
    guoduin = np.random.randint(0, guodu)
    chosen_RC[index[guoduin]] = 0
    guodu = code[index[guoduin]] - 0.5
    if guodu < 0:
        code[index[guoduin]] = 0
    else:
        code[index[guoduin]] = guodu

X[0, :RCc] = code

# 选择再制造商
chosen_MM = X[0, RCc + RCc * WPc + 1: RCc + RCc * WPc + MMc + 1]
code = chosen_MM.copy()
chosen_MM[chosen_MM > 0.5] = 1
chosen_MM[chosen_MM == 0.5] = 1
chosen_MM[chosen_MM < 0.5] = 0

while np.sum(chosen_MM) < q3:
    index = np.where(chosen_MM == 0)[0]
    guodu = len(index)
    guoduin = np.random.randint(0, guodu)
    chosen_MM[index[guoduin]] = 1
    guodu = code[index[guoduin]] + 0.5
    if guodu > 1:
        code[index[guoduin]] = 1
    else:
        code[index[guoduin]] = guodu

while np.sum(chosen_MM) > q4:
    index = np.where(chosen_MM == 1)[0]
    guodu = len(index)
    guoduin = np.random.randint(0, guodu)
    chosen_MM[index[guoduin]] = 0
    guodu = code[index[guoduin]] - 0.5
    if guodu < 0:
        code[index[guoduin]] = 0
    else:
        code[index[guoduin]] = guodu

X[0, RCc + RCc * WPc + 1: RCc + RCc * WPc + MMc + 1] = code

# 修正代码，避免工地向回收中心运输废料的量为零
guodu = X[0, RCc + 1: RCc + RCc * WPc + 1]
while True:
    guoduvalue = np.random.rand()
    if guoduvalue > 0:
        guodu[guodu == 0] = np.random.rand()
        break

X[0, RCc + 1: RCc + RCc * WPc + 1] = guodu

# 修正代码，避免回收中心向再生厂商运输废料的量为零
guodu = X[0, RCc + RCc * WPc + MMc + 1: RCc + RCc * WPc + MMc + MMc * RCc + 1]
while True:
    guoduvalue = np.random.rand()
    if guoduvalue > 0:
        guodu[guodu == 0] = np.random.rand()
        break

X[0, RCc + RCc * WPc + MMc + 1: RCc + RCc * WPc + MMc + MMc * RCc + 1] = guodu

X1 = X.copy()

# Uij 的废料运输量
guodu = X[0, RCc + 1: RCc + RCc * WPc + 1]
allocate_RC = np.reshape(guodu, (WPc, RCc))
rate = np.zeros((WPc, RCc))

for i in range(WPc):
    guodu = allocate_RC[i, :] * chosen_RC
    rate[i, :] = guodu / np.sum(guodu)

allocate_RC = rate
Uij = np.zeros((WPc, RCc))

for i in range(WPc):
    for j in range(RCc):
        Uij[i, j] = WP[i, 2] * allocate_RC[i, j]

# Ujk 的废料运输量
guodu = X[0, RCc + RCc * WPc + MMc + 1: RCc + RCc * WPc + MMc + MMc * RCc + 1]
allocate_MM = np.reshape(guodu, (RCc, MMc))
rate = np.zeros((RCc, MMc))

for j in range(RCc):
    guodu = allocate_MM[j, :] * chosen_MM
    rate[j, :] = guodu / np.sum(guodu)

allocate_MM = rate
Ujk = np.zeros((RCc, MMc))

for j in range(RCc):
    for k in range(MMc):
        Ujk[j, k] = epsilon_MM * RC[j, 6] * np.sum(Uij[:, j]) * allocate_MM[j, k]

# Ujm 的废料运输量
Ujm = np.zeros(RCc)

for j in range(RCc):
    Ujm[j] = epsilon_SM * RC[j, 6] * np.sum(Uij[:, j])

# Ujn 的废料运输量
Ujn = np.zeros(RCc)

for j in range(RCc):
    Ujn[j] = (1 - RC[j, 6]) * np.sum(Uij[:, j])

# Ukm 的废料运输量
Ukm = np.zeros(MMc)

for k in range(MMc):
    Ukm[k] = MM[k, 6] * np.sum(Ujk[:, k])

# Ukn 的废料运输量
Ukn = np.zeros(MMc)

for k in range(MMc):
    Ukn[k] = (1 - MM[k, 6]) * np.sum(Ujk[:, k])

# 计算直接碳排放量
def DIST(x1, y1, x2, y2):
    return np.sqrt((x1 - x2)**2 + (y1 - y2)**2)

TranCij = np.zeros((WPc, RCc))

for i in range(WPc):
    for j in range(RCc):
        TranCij[i, j] = Uij[i, j] * DIST(WP[i, 0], WP[i, 1], RC[j, 0], RC[j, 1]) * beta

TranCjk = np.zeros((RCc, MMc))

for j in range(RCc):
    for k in range(MMc):
        TranCjk[j, k] = Ujk[j, k] * DIST(RC[j, 0], RC[j, 1], MM[k, 0], MM[k, 1]) * beta

TranCjm = np.zeros(RCc)

for j in range(RCc):
    TranCjm[j] = Ujm[j] * DIST(RC[j, 0], RC[j, 1], SM[0], SM[1]) * beta

TranCjn = np.zeros(RCc)

for j in range(RCc):
    TranCjn[j] = Ujn[j] * DIST(RC[j, 0], RC[j, 1], LF[0], LF[1]) * beta

TranCkm = np.zeros(MMc)

for k in range(MMc):
    TranCkm[k] = Ukm[k] * DIST(MM[k, 0], MM[k, 1], SM[0], SM[1]) * beta

TranCkn = np.zeros(MMc)

for k in range(MMc):
    TranCkn[k] = Ukn[k] * DIST(MM[k, 0], MM[k, 1], LF[0], LF[1]) * beta

DCE = np.sum(TranCij) + np.sum(TranCjk) + np.sum(TranCjm) + np.sum(TranCjn) + np.sum(TranCkm) + np.sum(TranCkn)

# 计算间接碳排放量
CCij = np.zeros((WPc, RCc))

for i in range(WPc):
    for j in range(RCc):
        CCij[i, j] = Uij[i, j] * RC[j, 5]

CCjk = np.zeros((RCc, MMc))

for j in range(RCc):
    for k in range(MMc):
        CCjk[j, k] = Ujk[j, k] * MM[k, 5]

CCjm = np.zeros(RCc)

for j in range(RCc):
    CCjm[j] = Ujm[j] * SM[3]

CCjn = np.zeros(RCc)

for j in range(RCc):
    CCjn[j] = Ujn[j] * LF[3]

ICE = np.sum(CCij) + np.sum(CCjk) + np.sum(CCjm) + np.sum(CCjn)

TotalCC = DCE + ICE

# 成本计算
FC = np.sum(chosen_RC * RC[:, 3]) + np.sum(chosen_MM * MM[:, 3])

CopCij = np.zeros((WPc, RCc))

for i in range(WPc):
    for j in range(RCc):
        CopCij[i, j] = Uij[i, j] * RC[j, 4]

CopCjk = np.zeros((RCc, MMc))

for j in range(RCc):
    for k in range(MMc):
        CopCjk[j, k] = Ujk[j, k] * MM[k, 4]

CopCjm = np.zeros(RCc)

for j in range(RCc):
    CopCjm[j] = Ujm[j] * SM[2]

CopCjn = np.zeros(RCc)

for j in range(RCc):
    CopCjn[j] = Ujn[j] * LF[2]

VC = np.sum(CopCij) + np.sum(CopCjk) + np.sum(CopCjm) + np.sum(CopCjn)

TCij = np.zeros((WPc, RCc))

for i in range(WPc):
    for j in range(RCc):
        TCij[i, j] = Uij[i, j] * DIST(WP[i, 0], WP[i, 1], RC[j, 0], RC[j, 1]) * alpha

TCjk = np.zeros((RCc, MMc))

for j in range(RCc):
    for k in range(MMc):
        TCjk[j, k] = Ujk[j, k] * DIST(RC[j, 0], RC[j, 1], MM[k, 0], MM[k, 1]) * alpha

TCjm = np.zeros(RCc)

for j in range(RCc):
    TCjm[j] = Ujm[j] * DIST(RC[j, 0], RC[j, 1], SM[0], SM[1]) * alpha

TCjn = np.zeros(RCc)

for j in range(RCc):
    TCjn[j] = Ujn[j] * DIST(RC[j, 0], RC[j, 1], LF[0], LF[1]) * alpha

TCkm = np.zeros(MMc)

for k in range(MMc):
    TCkm[k] = Ukm[k] * DIST(MM[k, 0], MM[k, 1], SM[0], SM[1]) * alpha

TCkn = np.zeros(MMc)

for k in range(MMc):
    TCkn[k] = Ukn[k] * DIST(MM[k, 0], MM[k, 1], LF[0], LF[1]) * alpha

TC = np.sum(TCij) + np.sum(TCjk) + np.sum(TCjm) + np.sum(TCjn) + np.sum(TCkm) + np.sum(TCkn)

CTC = T * TotalCC

zeta = np.zeros(MMc)

for j in range(RCc):
    if RC[j, 6] > R2:
        zeta[j] = a2 * (R2 - RC[j, 6]) * np.sum(Uij[:, j])
    elif RC[j, 6] < R1:
        zeta[j] = a1 * (RC[j, 6] - R1) * np.sum(Uij[:, j])
    else:
        zeta