Fast Linearized Models

[Javier.Gil]

Hi all,

First of all I would like to thank you for such a great and well documented project on Wind Turbine simulation.

I was doing a bit of research on Wind Turbine dynamics and control. Unfortunately, I´m not an expert on control issues, so apologize in advance for the obvious questions on that field.

Taking the model NRELOffshrBsline5MW_Onshore, I have performed a linearization at several operating points in region 3 (above rated Wind Speed).

A standard header of one of my linearized solutions:

Order of States in Linearized State Matrices:

Row/column 1 = 1st tower fore-aft bending mode DOF (internal DOF index = DOF_TFA1)
Row/column 2 = 1st tower side-to-side bending mode DOF (internal DOF index = DOF_TSS1)
Row/column 3 = Variable speed generator DOF (internal DOF index = DOF_GeAz)
Row/column 4 = Drivetrain rotational-flexibility DOF (internal DOF index = DOF_DrTr)
Row/column 5 = 1st flapwise bending-mode DOF of blade 1 (internal DOF index = DOF_BF(1,1))
Row/column 6 = 1st flapwise bending-mode DOF of blade 2 (internal DOF index = DOF_BF(2,1))
Row/column 7 = 1st flapwise bending-mode DOF of blade 3 (internal DOF index = DOF_BF(3,1))
Row/column 8 = 1st edgewise bending-mode DOF of blade 1 (internal DOF index = DOF_BE(1,1))
Row/column 9 = 1st edgewise bending-mode DOF of blade 2 (internal DOF index = DOF_BE(2,1))
Row/column 10 = 1st edgewise bending-mode DOF of blade 3 (internal DOF index = DOF_BE(3,1))

Order of Control Inputs in Linearized State Matrices:

Column 1 = electrical generator torque (N·m) 4.30930E+04 op
Column 2 = rotor collective blade pitch (rad) 1.98206E-01 op
Column 3 = individual pitch of blade 1 (rad) 1.98206E-01 op
Column 4 = individual pitch of blade 2 (rad) 1.98206E-01 op
Column 5 = individual pitch of blade 3 (rad) 1.98206E-01 op

Order of Input Wind Disturbances in Linearized State Matrices:

Column 1 = horizontal hub-height wind speed (m/s) 1.60000E+01 op
Column 2 = horizontal wind direction (rad) 0.00000E+00 op
Column 3 = vertical wind speed (m/s) 3.47900E+00 op
Column 4 = horizontal shear parameter (-) 0.00000E+00 op
Column 5 = vertical power law shear exponent (-) 2.00000E-01 op
Column 6 = linear vertical shear parameter (-) 0.00000E+00 op
Column 7 = horizontal hub-height wind gust (m/s) 0.00000E+00 op

Order of Output Measurements in Linearized State Matrices:

Row 1 = GenSpeed (rpm)
Row 2 = GenTq (kN·m)
Row 3 = GenPwr (kW)
Row 4 = RotSpeed (rpm)
Row 5 = RotThrust (kN)
Row 6 = RotTorq (kN·m)
Row 7 = RotPwr (kW)
Row 8 = RootFxc1 (kN)
Row 9 = RootFxc2 (kN)
Row 10 = RootFxc3 (kN)
Row 11 = RootFyc1 (kN)
Row 12 = RootFyc2 (kN)
Row 13 = RootFyc3 (kN)
Row 14 = RootFzc1 (kN)
Row 15 = RootFzc2 (kN)
Row 16 = RootFzc3 (kN)
Row 17 = RootMxc1 (kN·m)
Row 18 = RootMxc2 (kN·m)
Row 19 = RootMxc3 (kN·m)
Row 20 = RootMyc1 (kN·m)
Row 21 = RootMyc2 (kN·m)
Row 22 = RootMyc3 (kN·m)
Row 23 = RootMzc1 (kN·m)
Row 24 = RootMzc2 (kN·m)
Row 25 = RootMzc3 (kN·m)
Row 26 = TipDxc1 (m)
Row 27 = TipDxc2 (m)
Row 28 = TipDxc3 (m)
Row 29 = TipDyc1 (m)
Row 30 = TipDyc2 (m)
Row 31 = TipDyc3 (m)
Row 32 = TipDzc1 (m)
Row 33 = TipDzc2 (m)
Row 34 = TipDzc3 (m)
Row 35 = TwrBsFxt (kN)
Row 36 = TwrBsFyt (kN)
Row 37 = TwrBsFzt (kN)
Row 38 = TwrBsMxt (kN·m)
Row 39 = TwrBsMyt (kN·m)
Row 40 = TwrBsMzt (kN·m)
Row 41 = TTDspFA (m)
Row 42 = TTDspSS (m)
Row 43 = TTDspAx (m)
Row 44 = BlPitch1 (deg)
Row 45 = BlPitch2 (deg)
Row 46 = BlPitch3 (deg)
Row 47 = HorWindV (m/sec)
Row 48 = HorWndDir (deg)
Row 49 = TSR (-)

In order to compare the performance of the linearized model with the nonlinear model, I prepared three basic simulations:

1.- Extreme Operating Wind Gust (max 8 m/s) in Closed Loop Condition.
2.- 2 Deg Pitch Increment in Open Loop Condition.
3.- 4 m/s Wind Speed Increment in Closed Loop Condition.

For the linearized model, I reproduce (as far as I can) the pitch control documented in:
Technical Report NREL/TP-500-38060: Definition of a 5-MW Reference Wind Turbine for Offshore System Development (great document, by the way).

In order to simplify the Variable Speed Control definition, I took the standard VS_CONT model with constant torque in region 3.

My first attempt was to use eigenanalysis.m routine and average out the linearized models obtained for different azimuts,
but then I found the so-clever MBC approach:
User’s Guide to MBC3 (Multi-blade Coordinate Transformation Utility for 3-Bladed Wind Turbines).

So, comparing the results of the FAST simultations with the MBC-Linearized model, following questions arise:

1.- Extreme Operating Wind Gust in Closed Loop Condition.

The linearized model exhibits a good agreement with FAST simulation, but I wonder if is there any way to reproduce the 1p-cyclic effects in the linearized model.
I realize that the cyclic effects are filtered out when performing the linearization, but maybe there exist an alternative way to deal with this effects. Otherwise,
I understand I could simulate the whole nonlinear model of the Wind Turbine in the Simulink environment with the aid of FAST_SFunc.mexw32.

It would be maravillous if one could develop individual pitch control strategies with the aid of the linearized models. Any Ideas?

2.- 2 Deg Pitch Increment in Open Loop Condition.

I have a problem with the linearized model. My thougth is that the problem is related with the evolution of the state variables.

What I expected to obtain was a steady condition in an operating point with a different rotational speed, but the offset error in the DOF_GeAz derivative causes non stationary terms in some output variables.
Mainly those which has a 1p periodic behaviour.

I try to work arround this problem, but no way to do that. I think it could maybe have to do with the procedure used to obtain the averaged matrices, which in this case is a simple average of all the matrices, let say A, B, Bd, C, D and Dd MBC transformed matrices.
From control point of view, this behaviour is related to the presence of integral terms, but I would need more understanding on the problem to filter out this non-desired terms.

This graph represents the evolution of RootFxc1 in this loadcase. As mentioned, it would be nice to obtain a steady value in this output.

See Figure PitchDelta_RootFxc1_(kN).png

3.- 4 m/s Wind Speed Increment in Closed Loop Condition.

In this case, the pitch controller command an increment of the collaborative pitch in order to recover the rated speed (also the rated power).

What I would spect in this case is a steady condition in all the outputs, but I don´t.

If we take a look at a variable as the Rotor Thrust, it´s OK.

See Figure WindSpeedDelta_RotThrust_(kN).png

If we analyze the evolution of RootFxc1, what we see is that a periodic behaviour arises.

See Figure WindSpeedDelta_RootFxc1_(kN).png

For this variable, the FAST model propose a steady situation with lower mean and higher amplitude so what I would spect for my LTI model was a constant steady value, beter than a periodic one.
On the other hand, it would be nice if the LTI model would be able to reproduce this kind of effects. I suppose this is a kind of MBC transformation of the offset error in one of the state variables.
Any idea about this topic?

Thank you in advance.

Kind Regards,

Javier.

[Jason.Jonkman]

Dear Javier,

Before answering your questions, I need to better understand what you’ve done. My understand is:
*You’ve used the linearization feature of FAST to obtain linear periodic models of the land-based NREL 5-MW turbine about several operating points.
*You’ve applied MBC3 to the linear periodic model, then azimuth-averaged the resulting matrices to obtain an LTI model.
*You’ve wrapped the LTI model inside a time-domain integrator in order to compare with responses obtained from FAST.
*For the closed-loop cases with the LTI model, you’ve mimicked the pitch controller of the NREL 5-MW in your time-domain integrator.
*For both the open- and closes-loop cases with the LTI model, you’ve replaced the torque controller of the NREL 5-MW with the simple variable-speed control option (VSContrl = 1) available in FAST.

Can you confirm my understanding of your approach?

Now here are my questions:
*Was the time integration of your LTI model and the control sytem implemented within MATLAB/Simulink?
*Have you added back in the periodic steady state values to the outputs (y vector) in order to compare to the outputs of FAST? (That is, y = yop + deltay.)
*Have you applied the inverse MBC transform to the outputs (y vector) in order to compare to the outputs of FAST? (Of course, this is only needed for the outputs effected by MBC, i.e. the blade motions and loads, such as RootFx1.)

Best regards,

[Javier.Gil]

Dear Jason,

I Apologize if I went too FAST in my first post. I´ll try to answer online to your questions:

*You’ve used the linearization feature of FAST to obtain linear periodic models of the land-based NREL 5-MW turbine about several operating points.

The model is correct: land-based NREL 5-MW turbine. As my aim was to control in region 3 of VS_CONT, I linearized the model for several horizontal hub-height wind speeds:
12:2:24 (m/s). In this case, plottings correspond to the 16 (m/s) model.

*You’ve applied MBC3 to the linear periodic model, then azimuth-averaged the resulting matrices to obtain an LTI model.

In this case, I used the mbc3 routines to obtain the MBC-averaged matrices. As mbc3 only computes the averaged state matrix, I added some lines to the script:

MBC_AvgB = MBC_AvgB + MBC_B(:,:,iaz);
if (exist('BMat'))
   MBC_AvgB = MBC_AvgB/NAzSteps;   % azimuth-average 
end

in order to obtain the rest of the averaged matrices. I wonder if there is a more efficient error-free averaging procedure for all these matrices.
I use:

2        MdlOrder

in my lin_options.dat file, but I think this point is not of relevance in this case.

*You’ve wrapped the LTI model inside a time-domain integrator in order to compare with responses obtained from FAST.

I prepared a simple wrapper in Simulink in order to integrate the MBC-averaged state-space representation of the nonlinear model. I used a rk4 fixed-step for integration.

*For the closed-loop cases with the LTI model, you’ve mimicked the pitch controller of the NREL 5-MW in your time-domain integrator.

As proposed in NREL/TP-500-38060, I used a gain scheduled PI controller. In the case of my last post, for 16 m/s wind speed an 11.36 Deg nominal pitch, I selected the corresponding scheduled gains.
I realize that some differences may arise due to gain changes in the nonlinear model, but I dont think this has something to do with the offset problems I documented earlier.

*For both the open- and closes-loop cases with the LTI model, you’ve replaced the torque controller of the NREL 5-MW with the simple variable-speed control option (VSContrl = 1) available in FAST.

For both LTI and FAST models, I used the same appoach: simple VS_CONT model:

   1        VSContrl   
 1173.0     VS_RtGnSp  
43093.0     VS_RtTq    
3.1318E-2   VS_Rgn2K   
9999.9E-9   VS_SlPc    

In the case of the simulink model, I wrote a embedded matlab function that emulates this behaviour. I realize that this model introduces a kind of non-linearity, and is maybe a part of the problem.
For VS_CONT controlling purposes, I use the derivative of DOF_GeAz, which I´m not sure is the best variable.

Can you confirm my understanding of your approach?

Now here are my questions:
*Was the time integration of your LTI model and the control sytem implemented within MATLAB/Simulink?

As said before, I build a simple simulink wrapper for the averaged MBC state-space representation of the model.

*Have you added back in the periodic steady state values to the outputs (y vector) in order to compare to the outputs of FAST? (That is, y = yop + deltay.)

No, I dont. For plotting purposes, I remove the operational values from the FAST nonlinear solution. My aim was to compare the performance of the models in terms of the solution mean values.
I know that I can obtain better approximations when adding operational values (y = yop + deltay) to the LTI solution, but I dont know how to do it properly when having rotor speed variations. As far as I understand, I could transform the 36 non-averaged yop values (each 10 Deg) into a look-up table and enter in this table with the azimut position. For MBC transformation purposes, I assume that azimut position can be approximated by:

azimut = (RotSpeed_op+deltaRotSpeed)*t

which I’m not sure is the best approach (see figure). Is this the correct way to do it?
It would be maravillous if it could be done in real-time. That way, one should be able to develop more efficient control schemes.

*Have you applied the inverse MBC transform to the outputs (y vector) in order to compare to the outputs of FAST? (Of course, this is only needed for the outputs effected by MBC, i.e. the blade motions and loads, such as RootFx1.)

For plotting purposes, I applied back the MBC transformation as follows:

function Ur = MBC_trans(Unr, rotpos)

% inputs:
% Unr		[3x1] non-rotating MBC set.
% rotpos	[1] azimut position. Calculated with the former formula.
% output:
% Ur		[3x1] rotating MBC set.

az1 = 0 + rotpos;
az2 = az1 + 2*pi/3;
az3 = az2 + 2*pi/3;

% compute transformation matrix

c1 = cos(az1);
c2 = cos(az2);
c3 = cos(az3);
s1 = sin(az1);
s2 = sin(az2);
s3 = sin(az3);

cos_col = [c1;c2;c3];
sin_col = [s1;s2;s3];

tt = [ones(3,1), cos_col, sin_col];

Ur = tt*Unr;

After checking RotTripletIndicesStates, RotTripletIndicesCntrlInpt, RotTripletIndicesOutput from mbc3.m, I assume that the MBC transformation is correct and the state space representation is in the non-rotating frame.
Therefore, I applied back the MBC transformation to the same RotTripletIndicesOutput.

Thank you in avance.

Kind regards,

Javier.

[Jason.Jonkman]

Dear Javier,

Thanks for clarifying your approach.

It is difficult to answer your questions because you have supplied a lot of information and are using more than just NREL-developed codes and models. Neither do I have the time to assess all of all your steps in detail. So, instead of answering your direct questions, I’ll try to give provide general information.

In general,the linear model should follow the nonlinear well for small perturbations about the linearization point. If your linear model is not tracking the nonlinear model well for small perturbations, I suspect there is a problem somewhere in the solution or comparison. It may help to simplify the model as much as possible in order to track down a problem (e.g., start by eliminating the generator azimuth DOF).

MBC is useful for rotors with 3 or more blades. MBC does not eliminate the periodicity, it simply transforms the rotating states of a model into the nonrotating frame (a linear model of mixed rotating and nonrotating states can then be represented in the nonrotating frame). For example, a 1p-cyclic load in the rotating frame contributes to 0p and 3p loading in the nonrotating frame. Much of the periodicity of the mixed frame linear matrices is captured by the MBC transform. Further azimuth-averaging of the lightly-periodic MBC-transormed matrices has been shown to be quite adequate, such that direct periodic techniques are likely not needed for 3-bladed rotors. See, for example, this reference: Stol A, Moll HG, Bir G, and Namik H. “A Comparison of Multi-Blade Coordinate Transformation and Direct Periodic Technicques for Wind Turbine Control Design.” 47th AIAA Aerospace Sciences Meeting Including The New Horizons Forum and Aerospace Exhibition, 5–8 January 2009, Orlando, FL. AIAA-2009-479. January 2009.

MBC is very useful to the design of independent blade pitch controllers. Again, see the reference above or review the papers it references.

Perhaps someone with more experience with the application of linear models and/or MBC can provide further guidance.

Best regards,

[Javier.Gil]

Dear Jason,

I will take a close look to the references you mentioned. Thanks a lot for your attention.

Perhaps you are wright when you say I must start with a simplifyed model, but before doing that, I have to make sure that I am not misunderstanding the FAST linearization process.

Let say I perform a linearization at an operation point that is above rated wind speed (pitch control zone). What I suppose is that the linearized state matrix refers to the dynamics of the open-loop wind turbine. Therefore, the contribution of selected inputs / disturbances are conveniently represented with [B] and [Bd] matrices.

Matrices [C], [D] and [Dd] allow us to build the open-loop transfer functions from selected inputs/disturbances to selected outputs, let say From BldPitch to RotSpeed. In this case, the only assumption I have to do is that matrices [B], [Bd], [C], [D] and [Dd] can be azimut-averaged in the same manner I buid the averaged state matrix.

In the case I analyze by example, 16 m/s wind speed, the results of the linearization for selected transfer functions (From GenTq, BldPitch and HorWindV to RotSpeed) are shown in the following figure:

My aim was to use all these averaged transfer functions in a simulink LTI model so as to start playing with different control algoritms. The point is that the system that I am trying to control, y = H(s) [u; ud],
appears to be somehow controlled, so I wonder if is there anything wrong with this basic approach.

Thank you in advance.

Kind Regards,

[Jason.Jonkman]

Dear Javier,

Yes, the linearized state matrix refers to the dynamics of the open-loop turbine (see exception below). And yes, the [C], [D] and [Dd] matrices refer to the open-loop dynamics from inputs/disturbances to outputs. And yes, the [B], [Bd], [C], [D] and [Dd] can be azimuth-averaged in the same manner that you build the averaged state matrix, [A]. Perhaps someone with more experience with linear systems analysis can comment on your Bode plot results.

Exception: the influence of the generator torque can be included in your linear model if the torque is implemented as a simple function of the turbine states (such as generator speed), without introducing new states. For example, the influence of the generator torque can be included when VSContrl = 1 or VSContrl = 0 with GenModel = 1 or 2. What torque option do you have set in your model before linearizing? Could your Bode plot results be influenced by the generator torque?

Best regards,

[Javier.Gil]

Dear Jason,

Thank you for your patience. As you suggested, I performed different simulations with different generator torque definitions (simple VSCONT and constant torque respectively), obtaining slightly different results. If I am wright, better results (also less damped) are obtained when using a constant torque definition, which is a normal assumtion when working in region 3. I also suppose that the differences are due to the non-linear behaviour of the simple VSCONT generator in regions 2 1/2 and 3.

The rest of the problems I documented earlyer had to do with numerical residuals in the state matrices. Therefore, I can now obtain the results that I expected. Moreover, I took my time in analyzing the evolution of the cyclic terms in the state matrices for different azimuth positions, so I get some insight on the effect of the MBC transformation. Thank you very much for the references you mentioned in your second post.

If I may be permitted to make one more (basic) question, I have one regarding to the tower shadow effect. If I am wright, the tower shadow has 3p contribution in non-rotating variables (ie. tilt moment in the yaw bearing) and 1p contribution in rotating variables (ie. Fx and My at the blade root of each blade). In order to evaluate this effect, I performed different simulations using the following parameters in the aerodyn input file:

90.000	HH
 0.3	TwrShad
 3.00	ShadHWid
 7.75	T_Shad_Refpt

(model: 5Mw baseline onshore turbine)

After analyzing results, I realized that the tower shadow effect is not so strong as I expected. Maybe the problem has to do with a bad selection of the former parameters. Maybe this tower shadow effect is somehow neglected when performing a linearization. Maybe I am misunderstanding the effect of the tower shadow, which I am afraid is the most probable thing.

The point is that I would like to make sure about which terms should be always taken in account and which terms can be negected.

Thank you in advance.

Kind regards,

[Jason.Jonkman]

Dear Javier,

I’m glad your MBC results now make sense to you.

AeroDyn v12.58 and earlier have only one tower-influence model. This is a tower-shadow model that only creates a wake deficit behind the tower, so it only influences the aerodynamics of a rotor operating downwind of the tower.

v13.00.00-bjj has two tower-influence models. The default model is the same tower-shadow model from v12.58, as described above. The second tower-influence model includes a tower dam model for upwind influence, as described in the AeroDyn Theory Manual: nrel.gov/docs/fy05osti/36881.pdf.

The inputs you’ve used apply only to the default tower-shadow model, and so, will not impact the aerodynamic loads of the NREL 5-MW baseline turbine, which has an upwind rotor.

Best regards,

[HamidReza.Esmaeilian]

Hello

As far as I know, FAST models the drivetrain of a wind turbine as follows:

2Ht d(ωt)/dt =Tw-Ktg(δt-δg)-Dtg(ωt-ωg)
d(δt-δg)/dt=d(δtg)/dt=ωt-ωg
2Hg d(ωg)/dt =Ktg(δt-δg)+Dtg(ωt-ωg)-Tg

where;
subscripts t and g indicate the rotor and generator, respectively;
ω: angular speed;
δtg: torsional angle;
Ktg: shaft stiffness;
Dtg: damping coefficient;
H: inertia constant;
Tg: electromagnetic torque;

thus we have three state variables, including Δωt, Δωg, and Δδtg.

In FAST linearization, I’d like to know which states model the mentioned variables.
Moreover, I’d like to know why GenDOF should be disabled when linearizing a constant-speed wind turbine. I think it leads to ignore variable Δωg.

regards

[Jason.Jonkman]

Dear Hamid Reza,

Your equations are effectively what FAST does for a 2-DOF generator + drivetrain model, except that I’m not sure where the "2"s in your equation come from.

For this 2-DOF model, FAST linearization will produce a 4-state model, including Δδg and Δωg corresponding to the generator DOF and Δδtg and Δωtg corresponding to the drivetrain-torsion DOF. As discussed in my Sep 02, 2014 response in the forum topic found here, it is common to remove the azimuth angle perturbation (Δδg) as a state. Then, with a little manipulation, you can transform this 4-state model into your desired 3-state model.

Best regards,

[Jason.Jonkman]

Dear Hamid Reza,

Regarding your second question, the generator DOF need not be disabled when linearizing a constant-speed wind turbine; in this case, it is possible to linearize such a model, keeping the slip of the generator.

Best regards,

[HamidReza.Esmaeilian]

Dear Jason

Thank you very much for your time and quick response

Just 1 more question:
you wrote that with a little manipulation, I can transform this 4-state model into my desired 3-state model. I was wondering if it means that I should remove the row and column corresponding to the azimuth angle perturbation (Δδg) from the state matrix (A).

thanks again.

[Jason.Jonkman]

Dear Hamid Reza,

Yes, removing the azimuth angle perturbation (Δδg) as a state involves removing a row and column from the state-space model. Then, with a little manipulation you can convert the resulting 3-state model from using Δωg, Δδtg, and Δωtg as states to using Δωt, Δωg, and Δδtg as states. That is: Δωt = Δωg + Δωtg.

Best regards,

[Joannes.Olondriz]

Dear All,

This is the most related post that I found it, but I still feeling confused.

Running GetMats.m function must I achieve wind turbine rotor blades dynamics expressed in rotating frames (AMat, BMat, …), while running MBC3.m function the dynamics of individual blades were expressed in a fixed (nonrotating) frame (MBC_A, MBC_B, …), right? And if I want to obtain more accurate model I should linearize the model at a number of different azimuth steps, right? Well, if I do that I need to take the average values to form my state space matrix (AvgAMat, AvgBMat,… in rotating reference frame), but:

• Why are not available this average value in non rotating reference frame (only I found MBC_AvgA)?
• Can I achieve the rest of this components (MBC_AvgB, MBC_AvgBd, MBC_AvgC, MBC_AvgD and MBC_AvgDd) as the same way that MBC_AvgA it is calculated?
  MBC_AvgA = zeros(ns);
  MBC_AvgA = MBC_AvgA + MBC_A(:,:,iaz);
  MBC_AvgA = MBC_AvgA/NAzSteps;

Thank you in advance.
Best regards,
JOE

[Jason.Jonkman]

Dear Joannes,

Yes, your understanding is correct. In general, I agree that AvgAMat etc. is less useful than MBC_AvgA etc.

Best regards,

[Joannes.Olondriz]

Dear Jason,

Thanks for yuor just in time reply.

So that, I would like to use MBC_Avg matrizes. But only is available MBC_AvgA. Thus to achieve the rest of matrizes I try the same as MBC_AvgA:

MBC_AvgB = sum(MBC_B(:,:,NAzimuths));
MBC_AvgB = MBC_AvgB/NAzSteps;

the same with Bd, C, D and Dd. But the result is not as expected:

In green: only 1 Azimuth and Mat
In reed: 4 Azimuths and AvgMat
In cyan: 4 Azimuths and MBC_Avg

I could say that the differences between the red and green is just the accuracy of applying 1 or 4 Azimuth steps, could I?
But why is such a big magnitude difference between red=green and cyan?
What can I do to achieve the average MBC matrix components (MBC_AvgB, MBC_AvgBd,…)?

Thank you so much.
Best regards,
JOE

[Joannes.Olondriz]

Dear Jason,

I fix it. Was my matrix managing error. Just when I did the division I didn’t take into account that I have to divide all vector elements.

So instead of doing:
MBC_AvgB = sum(MBC_B(:,:,NAzimuths));
MBC_AvgB = MBC_AvgB/NAzSteps;

It has to be done:
MBC_AvgB = sum(MBC_B(:,:,NAzimuths));
MBC_AvgB = MBC_AvgB./NAzSteps;

And the same for the rest Bd, C, D and Dd.

Best regards,
JOE