你如何为数组 w.r.t 标量声明部分，反之亦然？
How do you declare partials for arrays w.r.t scalars, and vice versa?
我正在尝试执行此计算：
使用以下代码：
class BendingMoment(om.ExplicitComponent):
def setup(self):
self.add_input('cableForce', units='N')
self.add_input('cableAngle', units='rad')
self.add_input('cablePosition', units='m')
self.add_input('FArray', units='N/m')
self.add_input('zArray', units='m')
self.add_input('wingspan', units='m')
self.add_output('MArray', units='N*m')
self.declare_partials('MArray', ['FArray', 'cableForce', 'cableAngle', 'cablePosition', 'wingspan', 'zArray'])
def compute(self, inputs, outputs):
Fc = inputs['cableForce']
ang = inputs['cableAngle']
FArray = inputs['FArray']
zArray = inputs['zArray']
b = inputs['wingspan']
a = inputs['cablePosition']
M = np.zeros_like(zArray)
for i in range(len(zArray)):
M[i] = np.trapz((zArray[i:]-zArray[i])*FArray[i:], zArray[i:]) - (0.5*b - a - zArray[i])*Fc*np.sin(ang)
outputs['MArray'] = M
其中Fc、ang、b、a是标量； FArray, zArray 是数组。
你如何为 M w.r.t 声明一个标量，即 Fc？
如果有一段代码获取一个数组并对它求和或积分以创建一个标量值，您将如何声明它的部分值？对于这两种情况，雅可比行列式究竟是什么样子的？
考虑大小为 (n0,n1) 的输出 Y 和大小为 (m0, m1) 的输入 X --- 两者都是二维数组，使事情更通用一些。然后，您向 OpenMDAO 提供的关于 Y 相对于 X 的导数的雅可比矩阵的大小将是一个大小为 (n0*n1, m0*m1) 的矩阵——即n0*n1 行和 m0*m1 列。
你得到的行数总是等于数组的展平大小（并且行的顺序与展平数组的顺序相同）。它对输入的工作方式类似。
查看偏导数雅可比行列式并查看您是否正确求导的最佳方法是使用 check_partials method，它将您的解析导数与近似导数进行比较。你可以使用它，即使你还没有定义你的导数——你只会得到所有的零，你可以通过在下面的代码中注释掉 compute_partials 来检查：
import numpy as np
import openmdao.api as om
class BendingMoment(om.ExplicitComponent):
def setup(self):
self.add_input('cableForce', units='N')
self.add_input('cableAngle', units='rad')
self.add_input('cablePosition', units='m')
self.add_input('wingspan', units='m')
self.add_input('FArray', units='N/m', shape=3)
self.add_input('zArray', units='m', shape=3)
self.add_output('MArray', units='N*m', shape=3)
self.declare_partials('MArray', ['FArray', 'cableForce', 'cableAngle', 'cablePosition', 'wingspan', 'zArray'])
def compute(self, inputs, outputs):
Fc = inputs['cableForce']
ang = inputs['cableAngle']
FArray = inputs['FArray']
zArray = inputs['zArray']
b = inputs['wingspan']
a = inputs['cablePosition']
M = outputs['MArray']
size = len(zArray)
# we can use np.trapz, only because we can "guess" how its implemented so that we can differentiate it
# for i in range(size):
# M[i] = np.trapz((zArray[i:]-zArray[i])*FArray[i:], zArray[i:]) - (0.5*b - a - zArray[i])*Fc*np.sin(ang)
# equivalent code to np.trapz written out with for loops explicitly. This code would be slower to run,
# so we'll leave the call to the numpy function, but we'll look at this to compute the partial derivatives
M[:] = 0 # clear out the array
for i in range(size):
for j in range(i, size-1):
M[i] += (zArray[j+1] - zArray[j]) * ((zArray[j] - zArray[i])*FArray[j] + (zArray[j+1] - zArray[i])*FArray[j+1])/2.
# # note that for i=j, you get a special case
# if i == j:
# M[i] += (zArray[j+1] - zArray[j]) * (zArray[j+1] - zArray[i])*FArray[j+1]/2.
# else:
# M[i] += (zArray[j+1] - zArray[j]) * ((zArray[j] - zArray[i])*FArray[j] + (zArray[j+1] - zArray[i])*FArray[j+1])/2.
M[i] -= (0.5*b - a - zArray[i])*Fc*np.sin(ang)
def compute_partials(self, inputs, J):
Fc = inputs['cableForce']
ang = inputs['cableAngle']
FArray = inputs['FArray']
zArray = inputs['zArray']
b = inputs['wingspan']
a = inputs['cablePosition']
size = len(zArray)
# We can't easily compute the partials for this component,
# because we don't know exactly what np.trapz does.
for i in range(size):
for j in range(i, size-1):
J['MArray', 'FArray'][i,j] += (zArray[j+1] - zArray[j]) * (zArray[j] - zArray[i])/2.
J['MArray', 'FArray'][i,j+1] += (zArray[j+1] - zArray[j]) * (zArray[j+1] - zArray[i])/2.
term1 = (zArray[j+1] - zArray[j])
term2 = ((zArray[j] - zArray[i])*FArray[j] + (zArray[j+1] - zArray[i])*FArray[j+1])/2.
# z portion
if i == j: # z == Array[j] -- special case
J['MArray', 'zArray'][i,i] += -term2 - term1 * (FArray[j+1])/2.
else:
J['MArray', 'zArray'][i,i] += -term1*(FArray[j+1]+FArray[j+1])/2.
J['MArray', 'zArray'][i,j] += - term2 + term1*FArray[j]/2.
J['MArray', 'zArray'][i,j+1] += term2 + term1*FArray[j+1]/2.
# account for the additional term with zArray in it
J['MArray', 'zArray'][i,i] += Fc*np.sin(ang)
J['MArray', 'cableAngle'][i] = -(0.5*b - a - zArray[i])*Fc*np.cos(ang)
J['MArray', 'cableForce'][i] = -(0.5*b - a - zArray[i])*np.sin(ang)
J['MArray', 'cablePosition'][i] = a*Fc*np.sin(ang)
J['MArray', 'wingspan'][i] = -0.5*b*Fc*np.sin(ang)
if __name__ == '__main__':
p = om.Problem()
p.model.add_subsystem('bending', BendingMoment(), promotes=['*'])
p.setup()
# note: leaving the zArray with the default value means it would be [1,1,1]
# which is non-physical, and will be a degenerate point to check the partials around.
# So we'll set it to a linear distribution
p['bending.zArray'] = np.linspace(0,10,3)
p.run_model()
p.check_partials()
给出此输出：
----------------------------------
Component: BendingMoment 'bending'
----------------------------------
bending: 'MArray' wrt 'FArray'
Forward Magnitude : 3.750000e+01
Fd Magnitude : 3.750000e+01 (fd:forward)
Absolute Error (Jfor - Jfd) : 6.343738e-09
Relative Error (Jfor - Jfd) : 1.691664e-10
Raw Forward Derivative (Jfor)
[[ 0. 25. 25. ]
[ 0. 0. 12.5]
[ 0. 0. 0. ]]
Raw FD Derivative (Jfd)
[[ 0. 24.99999999 25. ]
[ 0. 0. 12.5 ]
[ 0. 0. 0. ]]
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
bending: 'MArray' wrt 'cableAngle'
Forward Magnitude : 6.410044e+00
Fd Magnitude : 6.410039e+00 (fd:forward)
Absolute Error (Jfor - Jfd) : 4.991817e-06 *
Relative Error (Jfor - Jfd) : 7.787499e-07
Raw Forward Derivative (Jfor)
[[0.27015115]
[2.97166268]
[5.67317421]]
Raw FD Derivative (Jfd)
[[0.27015094]
[2.97166037]
[5.67316979]]
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
bending: 'MArray' wrt 'cableForce'
Forward Magnitude : 9.983052e+00
Fd Magnitude : 9.983052e+00 (fd:forward)
Absolute Error (Jfor - Jfd) : 3.700281e-09
Relative Error (Jfor - Jfd) : 3.706563e-10
Raw Forward Derivative (Jfor)
[[0.42073549]
[4.62809042]
[8.83544534]]
Raw FD Derivative (Jfd)
[[0.42073549]
[4.62809042]
[8.83544534]]
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
bending: 'MArray' wrt 'cablePosition'
Forward Magnitude : 1.457470e+00
Fd Magnitude : 1.457470e+00 (fd:forward)
Absolute Error (Jfor - Jfd) : 9.138199e-10
Relative Error (Jfor - Jfd) : 6.269904e-10
Raw Forward Derivative (Jfor)
[[0.84147098]
[0.84147098]
[0.84147098]]
Raw FD Derivative (Jfd)
[[0.84147099]
[0.84147099]
[0.84147099]]
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
bending: 'MArray' wrt 'wingspan'
Forward Magnitude : 7.287352e-01
Fd Magnitude : 7.287353e-01 (fd:forward)
Absolute Error (Jfor - Jfd) : 3.834701e-09
Relative Error (Jfor - Jfd) : 5.262132e-09
Raw Forward Derivative (Jfor)
[[-0.42073549]
[-0.42073549]
[-0.42073549]]
Raw FD Derivative (Jfd)
[[-0.4207355 ]
[-0.42073549]
[-0.42073549]]
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
bending: 'MArray' wrt 'zArray'
Forward Magnitude : 1.506254e+01
Fd Magnitude : 1.506254e+01 (fd:forward)
Absolute Error (Jfor - Jfd) : 9.989117e-07
Relative Error (Jfor - Jfd) : 6.631761e-08
Raw Forward Derivative (Jfor)
[[-9.15852902 0. 10. ]
[ 0. -4.15852902 5. ]
[ 0. 0. 0.84147098]]
Raw FD Derivative (Jfd)
[[-9.15852851 0. 10.00000049]
[ 0. -4.15852852 5.0000005 ]
[ 0. 0. 0.84147099]]
注意这里给出的代码是正确的，但不一定很快。一旦你有了正确的答案，向量化它们就容易多了。我发现首先使用速度更快的矢量化代码非常困难。我几乎总是将缓慢的 python-for-loop 版本的代码与矢量化代码
一起注释掉
我正在尝试执行此计算：
使用以下代码：
class BendingMoment(om.ExplicitComponent):
def setup(self):
self.add_input('cableForce', units='N')
self.add_input('cableAngle', units='rad')
self.add_input('cablePosition', units='m')
self.add_input('FArray', units='N/m')
self.add_input('zArray', units='m')
self.add_input('wingspan', units='m')
self.add_output('MArray', units='N*m')
self.declare_partials('MArray', ['FArray', 'cableForce', 'cableAngle', 'cablePosition', 'wingspan', 'zArray'])
def compute(self, inputs, outputs):
Fc = inputs['cableForce']
ang = inputs['cableAngle']
FArray = inputs['FArray']
zArray = inputs['zArray']
b = inputs['wingspan']
a = inputs['cablePosition']
M = np.zeros_like(zArray)
for i in range(len(zArray)):
M[i] = np.trapz((zArray[i:]-zArray[i])*FArray[i:], zArray[i:]) - (0.5*b - a - zArray[i])*Fc*np.sin(ang)
outputs['MArray'] = M
其中Fc、ang、b、a是标量； FArray, zArray 是数组。
你如何为 M w.r.t 声明一个标量，即 Fc？
如果有一段代码获取一个数组并对它求和或积分以创建一个标量值，您将如何声明它的部分值？对于这两种情况，雅可比行列式究竟是什么样子的？
考虑大小为 (n0,n1) 的输出 Y 和大小为 (m0, m1) 的输入 X --- 两者都是二维数组，使事情更通用一些。然后，您向 OpenMDAO 提供的关于 Y 相对于 X 的导数的雅可比矩阵的大小将是一个大小为 (n0*n1, m0*m1) 的矩阵——即n0*n1 行和 m0*m1 列。
你得到的行数总是等于数组的展平大小（并且行的顺序与展平数组的顺序相同）。它对输入的工作方式类似。
查看偏导数雅可比行列式并查看您是否正确求导的最佳方法是使用 check_partials method，它将您的解析导数与近似导数进行比较。你可以使用它，即使你还没有定义你的导数——你只会得到所有的零，你可以通过在下面的代码中注释掉 compute_partials 来检查：
import numpy as np
import openmdao.api as om
class BendingMoment(om.ExplicitComponent):
def setup(self):
self.add_input('cableForce', units='N')
self.add_input('cableAngle', units='rad')
self.add_input('cablePosition', units='m')
self.add_input('wingspan', units='m')
self.add_input('FArray', units='N/m', shape=3)
self.add_input('zArray', units='m', shape=3)
self.add_output('MArray', units='N*m', shape=3)
self.declare_partials('MArray', ['FArray', 'cableForce', 'cableAngle', 'cablePosition', 'wingspan', 'zArray'])
def compute(self, inputs, outputs):
Fc = inputs['cableForce']
ang = inputs['cableAngle']
FArray = inputs['FArray']
zArray = inputs['zArray']
b = inputs['wingspan']
a = inputs['cablePosition']
M = outputs['MArray']
size = len(zArray)
# we can use np.trapz, only because we can "guess" how its implemented so that we can differentiate it
# for i in range(size):
# M[i] = np.trapz((zArray[i:]-zArray[i])*FArray[i:], zArray[i:]) - (0.5*b - a - zArray[i])*Fc*np.sin(ang)
# equivalent code to np.trapz written out with for loops explicitly. This code would be slower to run,
# so we'll leave the call to the numpy function, but we'll look at this to compute the partial derivatives
M[:] = 0 # clear out the array
for i in range(size):
for j in range(i, size-1):
M[i] += (zArray[j+1] - zArray[j]) * ((zArray[j] - zArray[i])*FArray[j] + (zArray[j+1] - zArray[i])*FArray[j+1])/2.
# # note that for i=j, you get a special case
# if i == j:
# M[i] += (zArray[j+1] - zArray[j]) * (zArray[j+1] - zArray[i])*FArray[j+1]/2.
# else:
# M[i] += (zArray[j+1] - zArray[j]) * ((zArray[j] - zArray[i])*FArray[j] + (zArray[j+1] - zArray[i])*FArray[j+1])/2.
M[i] -= (0.5*b - a - zArray[i])*Fc*np.sin(ang)
def compute_partials(self, inputs, J):
Fc = inputs['cableForce']
ang = inputs['cableAngle']
FArray = inputs['FArray']
zArray = inputs['zArray']
b = inputs['wingspan']
a = inputs['cablePosition']
size = len(zArray)
# We can't easily compute the partials for this component,
# because we don't know exactly what np.trapz does.
for i in range(size):
for j in range(i, size-1):
J['MArray', 'FArray'][i,j] += (zArray[j+1] - zArray[j]) * (zArray[j] - zArray[i])/2.
J['MArray', 'FArray'][i,j+1] += (zArray[j+1] - zArray[j]) * (zArray[j+1] - zArray[i])/2.
term1 = (zArray[j+1] - zArray[j])
term2 = ((zArray[j] - zArray[i])*FArray[j] + (zArray[j+1] - zArray[i])*FArray[j+1])/2.
# z portion
if i == j: # z == Array[j] -- special case
J['MArray', 'zArray'][i,i] += -term2 - term1 * (FArray[j+1])/2.
else:
J['MArray', 'zArray'][i,i] += -term1*(FArray[j+1]+FArray[j+1])/2.
J['MArray', 'zArray'][i,j] += - term2 + term1*FArray[j]/2.
J['MArray', 'zArray'][i,j+1] += term2 + term1*FArray[j+1]/2.
# account for the additional term with zArray in it
J['MArray', 'zArray'][i,i] += Fc*np.sin(ang)
J['MArray', 'cableAngle'][i] = -(0.5*b - a - zArray[i])*Fc*np.cos(ang)
J['MArray', 'cableForce'][i] = -(0.5*b - a - zArray[i])*np.sin(ang)
J['MArray', 'cablePosition'][i] = a*Fc*np.sin(ang)
J['MArray', 'wingspan'][i] = -0.5*b*Fc*np.sin(ang)
if __name__ == '__main__':
p = om.Problem()
p.model.add_subsystem('bending', BendingMoment(), promotes=['*'])
p.setup()
# note: leaving the zArray with the default value means it would be [1,1,1]
# which is non-physical, and will be a degenerate point to check the partials around.
# So we'll set it to a linear distribution
p['bending.zArray'] = np.linspace(0,10,3)
p.run_model()
p.check_partials()
给出此输出：
----------------------------------
Component: BendingMoment 'bending'
----------------------------------
bending: 'MArray' wrt 'FArray'
Forward Magnitude : 3.750000e+01
Fd Magnitude : 3.750000e+01 (fd:forward)
Absolute Error (Jfor - Jfd) : 6.343738e-09
Relative Error (Jfor - Jfd) : 1.691664e-10
Raw Forward Derivative (Jfor)
[[ 0. 25. 25. ]
[ 0. 0. 12.5]
[ 0. 0. 0. ]]
Raw FD Derivative (Jfd)
[[ 0. 24.99999999 25. ]
[ 0. 0. 12.5 ]
[ 0. 0. 0. ]]
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
bending: 'MArray' wrt 'cableAngle'
Forward Magnitude : 6.410044e+00
Fd Magnitude : 6.410039e+00 (fd:forward)
Absolute Error (Jfor - Jfd) : 4.991817e-06 *
Relative Error (Jfor - Jfd) : 7.787499e-07
Raw Forward Derivative (Jfor)
[[0.27015115]
[2.97166268]
[5.67317421]]
Raw FD Derivative (Jfd)
[[0.27015094]
[2.97166037]
[5.67316979]]
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
bending: 'MArray' wrt 'cableForce'
Forward Magnitude : 9.983052e+00
Fd Magnitude : 9.983052e+00 (fd:forward)
Absolute Error (Jfor - Jfd) : 3.700281e-09
Relative Error (Jfor - Jfd) : 3.706563e-10
Raw Forward Derivative (Jfor)
[[0.42073549]
[4.62809042]
[8.83544534]]
Raw FD Derivative (Jfd)
[[0.42073549]
[4.62809042]
[8.83544534]]
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
bending: 'MArray' wrt 'cablePosition'
Forward Magnitude : 1.457470e+00
Fd Magnitude : 1.457470e+00 (fd:forward)
Absolute Error (Jfor - Jfd) : 9.138199e-10
Relative Error (Jfor - Jfd) : 6.269904e-10
Raw Forward Derivative (Jfor)
[[0.84147098]
[0.84147098]
[0.84147098]]
Raw FD Derivative (Jfd)
[[0.84147099]
[0.84147099]
[0.84147099]]
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
bending: 'MArray' wrt 'wingspan'
Forward Magnitude : 7.287352e-01
Fd Magnitude : 7.287353e-01 (fd:forward)
Absolute Error (Jfor - Jfd) : 3.834701e-09
Relative Error (Jfor - Jfd) : 5.262132e-09
Raw Forward Derivative (Jfor)
[[-0.42073549]
[-0.42073549]
[-0.42073549]]
Raw FD Derivative (Jfd)
[[-0.4207355 ]
[-0.42073549]
[-0.42073549]]
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
bending: 'MArray' wrt 'zArray'
Forward Magnitude : 1.506254e+01
Fd Magnitude : 1.506254e+01 (fd:forward)
Absolute Error (Jfor - Jfd) : 9.989117e-07
Relative Error (Jfor - Jfd) : 6.631761e-08
Raw Forward Derivative (Jfor)
[[-9.15852902 0. 10. ]
[ 0. -4.15852902 5. ]
[ 0. 0. 0.84147098]]
Raw FD Derivative (Jfd)
[[-9.15852851 0. 10.00000049]
[ 0. -4.15852852 5.0000005 ]
[ 0. 0. 0.84147099]]
注意这里给出的代码是正确的，但不一定很快。一旦你有了正确的答案，向量化它们就容易多了。我发现首先使用速度更快的矢量化代码非常困难。我几乎总是将缓慢的 python-for-loop 版本的代码与矢量化代码
一起注释掉