Several questions about OpenFAST Linearization

Dear all,

I encountered some difficulties during using OpenFAST Linearization. The version that I used is OpenFAST 3.4.1.

1. About TrimGain. When I performed Linearization analysis for OC4 DeepCWind offshore wind turbines, I found TrimGainoften lead to the oscillation in blade pitch angle. The TrimCase was set to 3. Hence, I read the related source code, and list them as follows:

select case (TrimCase_torque, TrimCase_pitch)
  xd%CtrlOffset = xd%CtrlOffset + (u%RotSpeed - p%RotSpeedRef) * p%TrimGain
end select
 if (p%TrimCase==TrimCase_pitch) then
      BlPitchCom = BlPitchCom + xd%CtrlOffset
end if

It seems that BlPitchCom is proportional to the integral of xd%CtrlOffset. xd%CtrlOffset is proportional to the integral of rotor speed error. Therefore, the transfer function from rotor speed error to blade pitch command contains Two integrators “1/s^2”. “1/s^2” in control field often results in bad dynamic response or system instability. I wonder whether it will be theoretically better to use only one integrator. For example, the source code might be changed into:

select case (TrimCase_torque, TrimCase_pitch)
  xd%CtrlOffset =  (u%RotSpeed - p%RotSpeedRef) * p%TrimGain  ! delete one xd%CtrlOffset  
end select
 if (p%TrimCase==TrimCase_pitch) then
      BlPitchCom = BlPitchCom + xd%CtrlOffset
end if

2.About TrimTol****. TrimTol is interpreted as the tolerance for the rotational speed convergence. I tried to find the detailed condition for convergence from the source code in order to better tune the value of TrimTol. Here is the related source code:

 !-①compute the error “eps_squared”
eps_squared = 0.0_ReKi
   do i = 1,p_FAST%Lin_NumMods 
      ThisModule = p_FAST%Lin_ModOrder( i ) !--indices that determine which order the modules are in the glue-code linearization matrix 
      do ThisInstance=1,size(y_FAST%Lin%Modules(ThisModule)%Instance)
        ny = y_FAST%Lin%Modules(ThisModule)%Instance(ThisInstance)%SizeLin(LIN_OUTPUT_COL) - y_FAST%Lin%Modules(ThisModule)%Instance(ThisInstance)%NumOutputs !last column before WriteOutputoccurs
       do j=1,ny
 indx = y_FAST%Lin%Modules(ThisModule)%Instance(ThisInstance)%LinStartIndx(LIN_OUTPUT_COL) + j - 1
             !LinStartIndx(LIN_OUTPUT_COL)
            if (EqualRealNos(m_FAST%Lin%y_interp( indx ), m_FAST%Lin%Y_prevRot( indx, m_FAST%Lin%AzimIndx ))) then !
               diff = 0.0_ReKi ! take care of some potential numerical issues
            else   
 diff = m_FAST%Lin%y_interp( indx ) - m_FAST%Lin%Y_prevRot( indx, m_FAST%Lin%AzimIndx )
            end if !
            eps_squared = eps_squared + ( diff / m_FAST%Lin%y_ref( indx ) ) ** 2         
end do!  
             end do    end do
   eps_squared = eps_squared / ( y_FAST%Lin%Glue%SizeLin(LIN_OUTPUT_COL) - y_FAST%Lin%Glue%NumOutputs )
  ②determine convergence.
 if (m_FAST%Lin%IsConverged) then!
         ! check that error equation is less than TrimTol !!!call 
         call calc_error(p_FAST, y_FAST, m_FAST, SrvD%y, eps_squared)
         m_FAST%Lin%IsConverged = eps_squared < p_FAST%TrimTol
      end if

“Convergence” at current azimuth is determined as follows. The differences for all LinOutputs at the same azimuth between this rotor revolution and the previous one is computed. y_ref (indx) is the difference between the maximum value and the minimum value for LinOutputs at the previous revolution. Then the eps_squared at current azimuth is defined as the mean value of the 2-norm of the differences for all LinOutputs at coterminous revolution with respect to the 2-norm of y_ref (indx). The linearization solution converges only if all eps_squared at NLinTimes azimuths are less than TrimTol. I wonder what is the difference between the output of “LIN_OUTPUT_COL” and the output of “y_FAST%Lin%Glue%NumOutputs”

3.About the perturbation of UserProp in Linearization. As shown in the source code, the perturbation of UserProp is 2pi/180. If the UserProp in airfoil files is measured in degree instead of rad and the UserProp is expected to change from -20 degree to 20 degree, I wonder whether it will be better to change p%du(28+k) from 2pi/180 to 2.

perturb = 2*D2R
 do k=1,p%NumBlades 
      p%du(28+k) = perturb                ! u%UserProp(:,:) = 29,30,31
 end do

4.As @jjonkman mentioned in issue #778 , enabling UserProp will enlarge the number of standard linearization inputs. I got the linearization results with UserProp. Elements from 10th to 47th in the linearization inputs are the AD User properties on blade 1, node 1-38. Elements from 48th to 85th are the AD User properties on blade 2, node 1-38. Elements from 86th to 123th are the AD User properties on blade 2, node 1-38. Now I would like to reduce the number of linearization inputs. I come up a method as follow: ①The first-order linearization state equation is : x = A xdot + Bu; y = Cx + Du; ②About the inputs: Assuming that the AD user properties at all nodes on one blade are the same, I consolidate the elements from 10th to 47th of inputs into one element, named “blade 1 UserProp”. Meanwhile, elements from 48th to 85th and elements from 86th to 123th are for “blade 2 UserProp” and “blade 3 UserProp”. Hence the inputs related to the AD UserProp are reduced from 38*3 to 3. ③About matrix B: For each row of B, I adds the elements from 10th to 47th and returns the result as the 10th element. The same handling is for the elements from 48th to 85th and the elements from 86th to 123th for each row of B. ④About matrix D: The handling of D is the same as B.

Hope to receive your kind response.

Yours

Liangwen Qi Chinese Academy of Science

Dear @ietqlw,

I'll answer regarding (1), (3), and (4), and let someone else comment on (2).

Regarding (1), I agree with you. The code matches the trim solution implementation plan (https://openfast.readthedocs.io/en/main/_downloads/34ccdeedadaf5bb6cc3af05369c33b8b/DevelopmentPlan-SetPoint-Linearization.pdf), but they are not correct for the intended proportional-only control, whereby the change in blade pitch should be proportional to the rotor speed error, not the integration of the error. This bug should be fixed.

Regarding (3), I agree as well. The internal perturbation size for UserProp assumes that UserProp is specified in radians, but if UserProp is specified in degrees, 2 may be a better perturbation size.

Regarding (4), I agree with how you plan to condense the nodel-level array of UserProp inputs into a single input per blade.

Best regards,

Regarding your question (2):

I wonder what is the difference between the output of “LIN_OUTPUT_COL” and the output of “y_FAST%Lin%Glue%NumOutputs”

The number of outputs in SizeLin(LIN_OUTPUT_COL) contains outputs for the OpenFAST input-output equations plus the number of user-requested outputs to write to the output file, NumOutputs. Since we don't want to include those user-requested outputs in the convergence criteria (and they are at the end of each module's list of outputs), we subtract them by module in the value for ny. Then, when we are dividing the eps_squared at the end, we take the total of all modules' outputs-- stored in Glue%SizeLin(LIN_OUTPUT_COL)--and subtract the total of all the modules' user-requested channels in the Glue%NumOutputs value. That denominator value, Glue%SizeLin(LIN_OUTPUT_COL) - Glue%NumOutputs could also be calculated by summing all of the ny values that are calculated above for each module, but that will give you the same answer as the way it is calculated in the code.

Dear @ietqlw,

I talked with @bjonkman about (1) and decided that we agree with you solution, but disagree with your statement that the original approach involving a double integration. The equation for BlPitchCom looks to be an integral; however, the BlPitchCom on the right-hand side of this expression is simply the pitch command defined by ServoDyn, which when linearization is enabled, is a fixed parameter (set to be the initial pitch, with no built-in integration). What we want is the pitch angle to be the initial pitch plus a perturbation proportional to the speed error. So, while we agree with you proposed solution, the result is that this change will eliminate the single integration that should not have been there in the first place.

Best regards,

Dear @jjonkman and @bjonkman,

Thank you for your clear answers and kindly sharing the trim solution implementation plan. There is no further doubt for questions(2)-(4).

Regarding question (1), I tested with my solution. The selected simulation object is the OC4 semi-submersible floating wind systems. The wind case is the steady wind type with the hub-height wind speed of 15 m/s and a wind shear coefficient of 0.14. However, I found that the convergence reached at the rotor speed larger than the reference one, as shown in the following screenshot.

blade pitch

rotor speed

I guess that the pure proportional trim of blade pitch angle leads to steady rotor speed error. If the initial blade pitch angle is set to zero, the pitch angle under the pure proportional trim is proportional to the rotor speed error. Then the blade pitch angle is programmed to zero as the speed tends to zero. But when blade pitch angle is zero, rotor speed will not tends zero. Hence, to get a steady operational point, a positive rotor speed error occurs to support a positive blade pitch angle.

Integral trim of blade pitch angle might help to eliminate the steady rotor speed error. Therefore, I change the source code to realize PI trim of blade pitch angle. Here are the changes:

(1)A variable DiscreteStateType%integer_spderr is added in ServoDyn_Registry.f90 to store the integer of the rotor speed error. typedef ^ DiscreteStateType ReKi integer_spderr - "integer of the rotor speed error qlw added" rad/s (2)Initialize DiscreteStateType%integer_spderr in SUBROUTINE SrvD_Init xd%integer_spderr = 0.0_ReKi !qlw added (3)Variables FAST_ParameterType%TrimIntGain, SrvD_InitInputType%TrimIntGain, SrvD_ParameterType%TrimIntGain are added in the Registry files. (4)Read FAST_ParameterType%TrimIntGain from the fast primary file in SUBROUTINE FAST_ReadPrimaryFile

CALL ReadVar( UnIn, InputFile, p%TrimIntGain    , "TrimIntGain", "qlw added: integral gain for the rotational speed error (>0) (rad/(rad/s) for yaw or pitch; Nm/(rad/s) for torque)", ErrStat2, ErrMsg2, UnEc)
      CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)

(5)Initialize SrvD_InitInputType%TrimIntGain in SUBROUTINE FAST_InitializeAll Init%InData_SrvD%TrimIntGain = p_FAST%TrimIntGain !qlw added (6)Initialize SrvD_ParameterType%TrimIntGain in SUBROUTINE SrvD_Init p%TrimIntGain = InitInp%TrimIntGain !qlw added (7)Add a channel for “TrimIntGain” in the OpenFAST primary input file. 0.0008965149 TrimIntGain - qlw added: integral gain for the rotational speed error (>0) (8)compute control offset xd%CtrlOffset in SUBROUTINE SrvD_UpdateDiscState.

real(Reki)                                     :: GKK_Trim     !qlw added
       select case (p%TrimCase)
case (TrimCase_torque, TrimCase_pitch)
     !    xd%CtrlOffset = xd%CtrlOffset + (u%RotSpeed - p%RotSpeedRef) * p%TrimGain
      GKK_Trim =  1.0/( 1.0 + u%BlPitch(1)/0.1099965 ) !gain-scheduling factor,0.1099965 rad = 6.302336 deg
          xd%integer_spderr = xd%integer_spderr + (u%RotSpeed - p%RotSpeedRef)*p%DT
          xd%CtrlOffset = (u%RotSpeed - p%RotSpeedRef) * p%TrimGain * GKK_Trim * 97.0_ReKi+ xd%integer_spderr * p%TrimIntGain * GKK_Trim * 97.0_ReKi  !97.0 is the gearbox ratio
      end select

Herein, I use the GS-PI presented in official document of the NREL 5MW baseline model. The proportional and integral gains at zero blade pitch angle are determined from the input channels “TrimGain” and “TrimIntGain” of primary input files. For OC4-semi wind turbines, I set them as:

0.006275604       TrimGain        - Proportional gain for the rotational speed error (>0) [used only if CalcSteady=True] (rad/(rad/s) for yaw or pitch; Nm/(rad/s) for torque)
0.0008965149    TrimIntGain     - qlw added: integral gain for the rotational speed error (>0) [used only if CalcSteady=True] (rad/(rad/s) for yaw or pitch; Nm/(rad/s) for torque)

Now, I have a few questions: Q1: Is there any potential problem in the above solution? Q2: Will the values of Twr_Kdmp and Bld_Kdmp affect the Linearization results? What are the reference setting values for them?

Dear @ietqlw,

Thinking about this more (I hadn't thought about the trim solution in quite awhile until you raised the question), it makes sense that the proportionally-only control results in a steady-state error of rotor speed for a second-order system. An integral term is needed to drive this error to zero. I was thinking the trim implementation was for a proportional-only control, but in reality it was implemented as an integral-only controller.

I haven't checked your source code changes, but I would expect a proportional-integral controller to convergence faster than an integral-only controller, as long as the proportional gain is set appropriately. But I would also say that the original code is OK as it is (but the documentation should refer to the gain as an integral gain rather than a proportional gain; this labeling of a proportional gain through me off).

Twr_Kdmp and Bld_Kdmp were briefly discussed in the following issue: https://github.com/OpenFAST/openfast/issues/853#issuecomment-940184369.

Best regards,

Dear @jjonkman,

Thank you for your reply. All of my doubts have been addressed.

Yours,

Liangwen Qi

OpenFAST / openfast

Several questions about OpenFAST Linearization #1545