扭转优化不起作用，因为初始向量已更改常量不保持常量

如何解决扭转优化不起作用，因为初始向量已更改常量不保持常量

我有以下代码框架，旨在优化由法线和双法线表示的扭转。这两个向量必须围绕轴旋转，切线。代码如下：

import numpy as np
from scipy.optimize import minimize_scalar
import vg
import math
from numpy import linalg as la
from pconst import const
PI = math.pi

def calc_torsion(iv,p,T,N,B):
    if (iv == p.shape[0]):
        Bp = B[1]
        Bz = B[0]
        pp = p[1]
        pz = p[0]
        Nz = N[0]
    elif (iv == p.shape[0]-1):
        Bp = B[0]
        Bz = B[iv]
        pp = p[0]
        pz = p[iv]
        Nz = N[iv]
    else:
        Bp = B[iv+1]
        Bz = B[iv]
        pp = p[iv+1]
        pz = p[iv]
        Nz = N[iv]
    Ba = Bp - Bz
    Bb = Bz - B[iv-1]
    a  = pp - pz
    b  = pz - p[iv-1]
    if (a[0]!=0 and b[0]!=0):
        BBx = 0.5*(Ba[0]/a[0] + Bb[0]/b[0]) 
    else:    
        BBx = 0
    if (a[1]!=0 and b[1]!=0):
        BBy = 0.5*(Ba[1]/a[1] + Bb[1]/b[1]) 
    else:    
        BBy = 0
    if (a[2]!=0 and b[2]!=0):
        BBz = 0.5*(Ba[2]/a[2] + Bb[2]/b[2])
    else:    
        BBz = 0    
    torsion = -(BBx*Nz[0] + BBy*Nz[1] + BBz*Nz[2])
    #print('Mittendrin: ',BBx,normal[0],BBy,normal[1],BBz,normal[2])
    return torsion

def rotate(v,a,s):
    nvec = np.dot(rotation_matrix(a,s),v)
    #print('Angle: ',angle,' Actual Angle: ',vg.angle(vector,nvec))
    if (np.absolute(s-np.deg2rad(vg.angle(v,nvec)))>0.001 and np.absolute(s+np.deg2rad(vg.angle(v,nvec))-2*PI)>0.001 and np.absolute(s+np.deg2rad(vg.angle(v,nvec)))>0.001):
        print('Angle differene in rotate(),Angle: ',s,np.deg2rad(vg.angle(v,nvec)))
    #print('Angle: ',theta,'norm Axis:',la.norm(axis),'norm Vec Pre:',prenormvec,'norm Vec Post:',la.norm(nvec))
    nvec /= la.norm(nvec) 
    return nvec    

def rotation_matrix(axis,theta):
    """
    Return the rotation matrix associated with counterclockwise rotation about
    the given axis by theta radians.
    """
    axis = np.asarray(axis)
    axis = axis / math.sqrt(np.dot(axis,axis))
    a = math.cos(theta / 2.0)
    b,c,d = -axis * math.sin(theta / 2.0)
    aa,bb,cc,dd = a * a,b * b,c * c,d * d
    bc,ad,ac,ab,bd,cd = b * c,a * d,a * c,a * b,b * d,c * d
    return np.array([[aa + bb - cc - dd,2 * (bc + ad),2 * (bd - ac)],[2 * (bc - ad),aa + cc - bb - dd,2 * (cd + ab)],[2 * (bd + ac),2 * (cd - ab),aa + dd - bb - cc]])
def opt_frenetframe(p,B):
    for iv in range(1,2):
        const.Tzero = T[iv]
        const.Nzero = N[iv]
        const.Bzero = B[iv]
        #Tzero = np.array([T[iv][0],T[iv][1],T[iv][2]])
        #Nzero = np.array([N[iv][0],N[iv][1],N[iv][2]])
        #Bzero = np.array([B[iv][0],B[iv][1],B[iv][2]])
        print('First  Tzero: ',const.Tzero,'  T[iv]: ',T[iv])
        print('First  Nzero: ',const.Nzero,N[iv])
        print('First  Bzero: ',const.Bzero,'  B[iv]: ',B[iv])
        print('------------------------------------------------------------------------------------------------')
        pre_torsion = np.absolute(calc_torsion(iv,B))
        #print('Point: ',iv,'Pre Torsion: ',pre_torsion)
        
        def f(theta):
            N[iv] = rotate(const.Nzero,theta)
            B[iv] = rotate(const.Bzero,theta)
            print('Tzero: ',T[iv])
            print('Nzero: ',N[iv])
            print('Bzero: ',B[iv])
            print('------------------------------------------------------------------------------------------------')
            t  = np.absolute(calc_torsion(iv,B))
            return t

        res = minimize_scalar(f)#,bounds=(-PI,PI),method='bounded')
        angle = res.x
        N[iv] = rotate(const.Nzero,angle)
        B[iv] = rotate(const.Bzero,angle)
        print('Last  Tzero: ',T[iv])
        print('Last  Nzero: ',N[iv])
        print('Last  Bzero: ',B[iv])
        post_torsion = np.absolute(calc_torsion(iv,B))
        print('Point: ',pre_torsion,'Angle: ','Post Torsion: ',post_torsion)
    return (T,B)

您可以使用以下几行代码开始：

p = np.array(((-0.48084623,-1.12269079,-0.10674437),(-0.48104423,-1.12605479,-0.05558529),(-0.48097343,-1.12661079,-0.00442211)))
T = np.array(((-0.01030932,-0.09440498,0.9954805),(-0.00124487,-0.03830994,0.99926513),(8.97378314e-05,1.50947006e-02,9.99886064e-01)))
N = np.array([[0.38339887,-0.91982235,-0.08325952],[0.38577626,-0.92193337,-0.03486459],[0.38697396,-0.92198611,0.01388396]])
B = np.array([[0.92352532,0.38080775,0.0456775],[0.92259153,0.38544936,0.01592676],[0.92209064,0.38692862,-0.00592399]])
T,B = opt_frenetframe(p,B)

问题是 Nzero 和 Bzero 在角度优化期间不会保持恒定，这就是优化不起作用的原因。当我使用

const.Tzero = np.array([T[iv][0],T[iv][2]])
const.Nzero = np.array([N[iv][0],N[iv][2]])
const.Bzero = np.array([B[iv][0],B[iv][2]])

代替

const.Tzero = T[iv]
const.Nzero = N[iv]
const.Bzero = B[iv]

代码运行良好，但为什么会这样呢？我用 id() 检查了内部地址，它们在 Nzero 和 N[iv]、Bzero 和 B[iv] 之间是不同的。那么，如果它不能是 Nzero 和 Bzero，因为它们是常数，那么相应地址处的值会发生什么变化？

我花了很长时间才找到代码中的错误，但我仍然不明白是什么导致了错误。预先非常感谢您。