Search code examples
openmdao

How do you declare partials for arrays w.r.t scalars, and vice versa?

I'm trying to implement this calculation:

enter image description here

With the following code:

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

Where Fc, ang, b, a are scalars; FArray, zArray are arrays.

How would you declare partials for M w.r.t a scalar, i.e Fc?

And if there was a code that took an array and summed or integrated it to create a scalar value, how would you declare the partials for that? What exactly would the Jacobian look like for both of these cases?

Solution

  • Consider an output Y of size (n0,n1) and an input X of size (m0, m1) --- both are 2D arrays to make things a bit more general. Then the size of the jacobian that you provide to OpenMDAO for derivatives of Y with respect to X will be a matrix of size (n0*n1, m0*m1) ---i.e. n0*n1 rows and m0*m1 columns.

    You always get the number of rows equal to the flattened size of the array (and the rows will be in the same order as the flattened array). It works similarly for inputs.

    The best way to look at the partial derivative Jacobian, and to see if you got the derivative right is to use the check_partials method which will compare your analytic derivatives to the approximate ones. You can use this, even if you haven't defined your derivatives yet --- you'll just get all zeros, which you can check by commenting out the compute_partials in the following code: 

    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()

    gives this output: 

    ----------------------------------
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]]

    Note that the code given here is correct, but not necessarily fast. Once you have correct answers, it is a lot easier to vectorize them. I find it very difficult to work with the faster, vectorized code first. I almost always keep the slow python-for-loop version of the code commented out along side the vectorized code