Feature request: Power flow visualization function

[dtgaebe]

Feature description. This function would take the optimization results and the WEC hydrodynamic properties and to illustrate the average power flows from wave through body to electrical power.

Describe the solution you'd like This function would quickly visualize where performance gains have been made, or where there is room for improvement.

For the 1 dof case, after Falnes 6.2 For a single body oscillating in one dof with complex velocity $U$, with intrinsic impedance $Z_i$, radiation resistance $R_i = Re(Zi)$ and excitation force coefficient $F{ex}$ the excitation power from the incident waves is given by (Falnes 6.3) $$P{ex} = \frac{1}{4} (F{ex}U^ + F_{ex}^U) $$ and the power radiated (away) by the bodies oscillation by (Falnes 6.4) $$P_{rad} = \frac{1}{2} Ri |U|^2$$ The maximal absorbable power is (Falnes 6.10) $$P{max} = \frac{|F_{ex}|^2}{8R_i}$$

The power absorbed by the WEC is then (Falnes 6.2) $$P{mech} = P{ex} -P{rad} $$ Which needs to equal the solver results for the mechanical power, thus we can use an assert statement. The power loss in the PTO can be calculated based on the solver results $$P{PTO,loss} = P{mech} -P{elec} $$

~~The power that has not been absorbed shall be the difference between the maximal absorbable power and the exited power $$P{not abs} \ne P{max} -P_{ex} $$~~

P_mech = pto.mechanical_average_power(wec, x_wec, x_opt, waves, nsubsteps)
P_elec = pto.average_power(wec, x_wec, x_opt, waves, nsubsteps)

hd = wot.linear_hydrodynamics(bem_data, mass, stiffness)
hd = wot.check_linear_damping(hd)
Zi = wot.hydrodynamic_impedance(hd)
Rad_res = np.real(Zi.squeeze())

Fex = wec_fdom.force.sel(type=['Froude_Krylov', 'diffraction']).sum('type')
Vel = wec_fdom.vel

P_excited_6_3 = (1/4)*(Fex*np.conjugate(Vel) + np.conjugate(Fex)*Vel ).squeeze().sum('omega').item()
P_radiated_6_4 = ((1/2)*(Rad_res * np.abs(Vel)**2 ).squeeze().sum('omega').item())
P_max_absorbable_6_10 = (np.abs(Fex)**2/(8*Rad_res) ).squeeze().sum('omega').item()
P_PTO_absorbed_6_2 = P_excited_6_3 - P_radiated_6_4
assert(np.isclose(P_PTO_absorbed_6_2, -1*P_mech)) #assert that solver solution matches
P_PTO_loss = P_mech - P_elec
P_not_absorbed = P_max_absorbable_6_10 - P_excited_6_3 

Power_flows = [P_max_absorbable_6_10 , P_PTO_loss, -1*P_radiated_6_4, -1*P_not_absorbed, P_elec, ]

As visualization of the different average power values could be done with a sankey diagram, for example like this

from matplotlib.sankey import Sankey
fig = plt.figure(figsize = [6,4])
ax = fig.add_subplot(1, 1, 1,)
sankey = Sankey(ax=ax, 
                scale=0.5/P_max_absorbable_6_10 ,
                offset= 0,
                format = '%d',shoulder = 0.02)

sankey.add(flows=Power_flows, 
           labels = ['Max Absorbable', 'PTO-Loss', 'Radiated', 'Non-absorbed', 'Elec'], 
           orientations=[0, -1, -1, -1, 0],#arrow directions,
           pathlengths = [0.4,0.15,0.2,0.1,0.3,],
           trunklength = 0.5,
           edgecolor = '#2b8cbe',
           facecolor = '#2b8cbe' )

diagrams = sankey.finish()
for diagram in diagrams:
    for text in diagram.texts:
        text.set_fontsize(10);
plt.axis("off") 
plt.show()

EDIT: add python syntax highlighting EDIT2: strike through wrong statement

[jtgrasb]

Currently, something is wrong with the implementation. In some cases, the excited power is shown as larger than the max absorbable power.

Regular wave: Amplitude = 0.125 m, Period = 2.5:

In this case, the non-absorbed power becomes negative, which shouldn't make sense. It seems to be either an issue with the equation for excited power or for max absorbable power. @ryancoe @cmichelenstrofer

[ryancoe]

@dtgaebe - I can't immediately find a flaw in your logic, but what @jtgrasb shows is obviously a problem.

@jleonqu @dforbush2 - Can you guys take a look at how @dtgaebe has proposed to build this power accounting Sankey diagram as you have both done similar things in the past (see, e.g., Fig. 12 in this paper: https://doi.org/10.1109/TSTE.2023.3237125)?

[dtgaebe]

I agree.

The maximal absorbable power (Falnes 6.10) requires the assumption of an optimal velocity $U{opt} = F{ex}/(2R_i) $. We have shown in the IEEE paper that the optimal velocity for optimizing the electrical power does not necessarily equal this optimal velocity (for mechanical power transfer). Instead, we require a WEC velocity that ensures that the generator current is in phase with the Thevenin voltage (Excitation of the WEC-PTO eequivalent system). The electrically optimal WEC velocity might require a non-negligible amount of reactive power (which might make the excitation power equation 6.3 invalid?)

I am wondering if we need to have another inflow - the WEC reactive power to not to violate conservation of energy?? @gbacelli

[dforbush2]

The last guys I did this for was CPower. I can provide the code for how we calculated all of these constituent terms for an idealized reactive-power controller...but it is a labyrinth of sub-functions so lmk if you'd like to subject yourself to that. The rub was in this case that

P_max_absorbable_6_10 = (np.abs(Fex)**2/(8*Rad_res) ).squeeze().sum('omega').item()

wasn't valid because their system was relied on the relative motion between two bodies, whereas this equation is only valid for a single body (or something like that...ask GB)

[dforbush2]

The other really tricky part here was making sure that electrical losses were calculated appropriately based on their definitions of R and Kt. Need to know if these are based off of RMS, phase-to-phase, phase-to-neutral, etc. I have a great summary from Steve on all this stuff.

I've attached the calculation from some time-domain simulink outputs, along with the relevant simulation. Main code is constrainedGains..., plotting code is makeEnergyFlowV2. In the main code see lines 236 through 253 and 553+. (Note the variable waveMax doesn't end up getting used). I had all the subfunctions at the end of the code because I apparently hated myself when I wrote this. WOTupload.zip

Checkout wavePowOut calculations, and note that viscous losses were calculated from the balance of the other power losses.

[dtgaebe]

I misunderstood this, but the max absorbed power actually can reach the PTO and shouldn't be all the way to the left of the plot.

Also in Falnes 6.10 he clearly states that the maximal absorbed power is half of the optimal incident excitation power and in this case equals the radiated power $$P{max} = \frac{|F{ex}|^2}{8Ri} =P{rad}^{opt} = \frac{1}{2}P{ex}^{opt} $$ So in order to maximize the mechanical power, equal amount of power needs to be radiated away. So my equation above calculating the not absorber power is definitely wrong! It should be $$P{ex}^{opt} - P_{ex}$$ We might have to reconsider how we illustrate the entire diagram. I think the quantities that we are actually calculating should be in an envelope of optimal power transfer, instead of trying to relate the quantities we're calculating with the optima.

The inside arrows would be our actual calculations, the outside arrows (green, purple, blue, light green and brown) the absolute optima. However, with this we have to be careful with communicating, since in practice it might not be desirable to entirely fill out the arrows... We can likely extend this on the electrical side as well, using the optimality conditions for the power transfer on the electrical side and then enclosing our optimization results

@jtgrasb for now an easy fix, that is definitely valid is to just focus on the actual results and not the upper bounds, but use the scale in the sankey diagram based on the max power, so that different cases are comparable:

[jleonqu]

@dtgaebe I agree with you, your first equation to calculate the not-absorbed power is not correct because the excitation power $P{ex}$ could be greater than the the maximum absorbable power $P{max}$, this is why @jtgrasb is having some negative results for the not-absorbed power. Your second equation, which is $P^{opt}{ex} - P{ex}$, should be equal to $P{max}-P{mech}$.

For the paper that @ryancoe mentioned above (available here), I used the following equations to calculate the power losses:

• Power loss in wave-body interaction: $P{max}-P{mech}$
• PTO losses: $P{mech} - P{elec}$

[dtgaebe]

Hmm, @jleonqu why do you think $P{ex}^{opt} - P{ex} $, should be equal to $P{max} - P{mech}$?

I think the "power loss" in wave body interactions should be the radiated power $$P_{rad} = \frac{1}{2} Ri |U|^2$$ the $P{max} - P_{mech}$ is to my understanding more an unused potential, I call it headroom in the more detailed graphic below.

[jtgrasb]

@dtgaebe Thanks for providing some diagrams - they've helped my understanding of the power flow significantly. For the EWTEC paper that I'm writing, I think it may actually be better to visualize the maximum absorbable power rather than the excited power. In this case, my input is the max absorbable power and I measure mechanical loss between the waves and WEC (P_loss,Wave-WEC = P_max - P_mech), the friction loss from the drivetrain friction (P_loss,fric = F_fric*vel), and the transfer losses between the WEC and PTO as well as PTO and generator (P_loss,transfer = P_mech - P_fric - P_elec). I think this representation better shows how the negative spring can improve the transfer of energy between the wave and WEC as well as improving the transfer of energy from WEC through the PTO even though friction losses are increased due to larger velocities. Do you think this makes sense? And do you have any suggestions for the naming convention?

[dtgaebe]

I like the idea of splitting up the losses!

But, to my current understanding calling $P{max} - P{mech}$ something with "loss" is wrong, since in $P{max} - P{mech}$ there can be plenty of power that has not been transferred to the WEC in the first place. The mechanical power that is actually lost on the hydrodynamic side is the power radiated away $P_{rad} = \frac{1}{2} R_i |U|^2$

Also, "transfer" is probably not the best term, since I assume that the losses closes to the electrical side are mostly due to the generator winding resistance (I just did a quick check and it's not super simple to make it a sum 1 game in WecOptTool, because we're not using the zero-frequency component in the pto impedance, so some losses are due to transfer, but most are due to I^2R)

[dtgaebe]

@jtgrasb

I extended the graph with winding loss and friction loss, which you can calculate with the real part of the drivetrain_impedance and winding_impedance $R{11}$ and $R{22}$, respectively.

Zi = wot.hydrodynamic_impedance(hd)
Rad_res = np.real(Zi.squeeze())

Fex = wec_fdom.force.sel(type=['Froude_Krylov', 'diffraction']).sum('type')
Vel = wec_fdom.vel

P_radiated_6_4 = ((1/2)*(Rad_res * np.abs(Vel)**2 ).squeeze().sum('omega').item())
P_excited_6_3 = (1/4)*(Fex*np.conjugate(Vel) + np.conjugate(Fex)*Vel ).squeeze().sum('omega').item()
P_PTO_absorbed_6_2 = P_excited_6_3 - P_radiated_6_4
assert(np.isclose(P_PTO_absorbed_6_2, -1*P_mech)) #assert that solver solution matches

#adding zero component
import xarray as xr
winding_res = xr.DataArray(np.real(pto_impedance_22)[0], dims=["omega"], coords=[("omega", np.concatenate(([0], hd.omega.values))),],)
drive_train_res = xr.DataArray(np.real(-1*pto_impedance_11)[0], dims=["omega"], coords=[("omega", np.concatenate(([0], hd.omega.values))),],)

P_loss_winding = (1/2)*(winding_res*np.abs(current)**2).squeeze().sum('omega').item()
P_loss_friction = (1/2)*(drive_train_res*np.abs(Vel)**2).squeeze().sum('omega').item()

assert(np.isclose(P_loss_friction+ P_loss_winding, -1*(P_mech-P_elec))) #assert that solver solution matches

Power_flows = [P_excited_6_3,  -1*P_loss_winding, -1*P_loss_friction, -1*P_radiated_6_4, P_elec, ]

fig = plt.figure(figsize = [6,4])
ax = fig.add_subplot(1, 1, 1,)
sankey = Sankey(ax=ax, 
                scale=0.5/P_max_absorbable,
                offset= 0,
                format = '%.2f',shoulder = 0.02)

sankey.add(flows=Power_flows, 
           labels = ['Exitation',  'PTO-elec-res', 'PTO-friction', 'Radiated', 'Elec'], 
           orientations=[0, -1, -1, -1, 0],#arrow directions,
           pathlengths = [0.4,0.1,0.2,0.2,0.3,],
           trunklength = 0.5,
           edgecolor = '#2b8cbe',
           facecolor = '#2b8cbe' )

diagrams = sankey.finish()
for diagram in diagrams:
    for text in diagram.texts:
        text.set_fontsize(10);
plt.axis("off") 
plt.show()

[jtgrasb]

Thanks @dtgaebe! I think it'll be very helpful to show both of those components!

[dforbush2]

I've been following this dialog and I want to caution you once again or maybe add a comment in the code that the expressions in the last piece of code are valid only for particular definitions of resistance, current, torque constant (if there is one...?), etc and it's worth noting the specific case this is for. These conversations get complicated and circular otherwise.

EDIT: deleted email attachement

[dtgaebe]

@dforbush2 thank you! Yes, this is for the WaveBot, with it's linear PTO impedance which @jtgrasb currently uses. The assumption behind the the 2-port PTO impedance is that it's a perfect gyrator between current and torque and vice versa (thus no power loss, thus the torque constant doesn't show up). More general the losses should still hold:

P_{loss}^{PTO} = Re\left(\begin{bmatrix}Z_{11} & Z_{12} \\ Z_{21} & Z_{22} \end{bmatrix} \right) \begin{bmatrix} |flow_1|^2 \\ |flow_2|^2 \end{bmatrix}

All the max values only exist in an unconstrained world... Careful with more degrees of freedom, or relative body motion. But Falnes' book does have theory about those cases as well. All quantities I'm describing are the Fourier series coefficients, but they do match the time series power outputs from the solver in the linear, unconstrained WaveBot case.

[dtgaebe]

Ok, @jtgrasb I just had the chance to talk to Giorgio. He recommended to keep it simple, but we need to emphasize that in order to have maximum absorption we need maximum radiation too. In most cases (suboptimal velocity) there will be a lot of unused potential. Giorgio previously called this reflected (following micro wave language), but we're concerned that this can be misunderstood in the ocean wave field.

People have to understand that max power transfer to the WEC requires 50% radiated power, which needs to go in the text, but I also added it to the figure (somewhat).

The biggest change that we should have right now for the EWTEC paper is that we start with the optimal excitation power at the very left! We can refined the composition of the Unused potential and the PTO-loss later.

P_max_absorbable_6_10 = (np.abs(Fex)**2/(8*Rad_res) ).squeeze().sum('omega').item()
P_opt_excitation_6_10 = 2*P_max_absorbable_6_10
P_radiated_6_4 = ((1/2)*(Rad_res * np.abs(Vel)**2 ).squeeze().sum('omega').item())
P_excited_6_3 = (1/4)*(Fex*np.conjugate(Vel) + np.conjugate(Fex)*Vel ).squeeze().sum('omega').item()
P_PTO_absorbed_6_2 = P_excited_6_3 - P_radiated_6_4
assert(np.isclose(P_PTO_absorbed_6_2, -1*P_mech)) #assert that solver solution matches
P_PTO_loss = P_mech - P_elec
P_unused = P_opt_excitation_6_10 - P_excited_6_3

Power_flows = [P_opt_excitation_6_10, P_PTO_loss, -1*P_radiated_6_4, -1*P_unused, P_elec, ]

fig = plt.figure(figsize = [6,4])
ax = fig.add_subplot(1, 1, 1,)
sankey = Sankey(ax=ax, 
                scale=0.5/P_max_absorbable,
                offset= 0,
                format = '%.2f W',shoulder = 0.02)

sankey.add(flows=Power_flows, 
           labels = ['Optimal Excitation \n $ 2 \\frac{|F_e|^2}{8Re(Z_i)} = 2*P_{mech}^{max}$ ', 
                      'PTO-Loss \n $ P_{mech} - P_{elec}$', 
                      'Radiated \n $ \\frac{1}{2} Re(Z_i) |U|^2  $ ', 
                      'Unused Potential \n(neither absorbed nor radiated)', 
                      'Electrical'], 
           orientations=[0, -1, -1, -1, 0],#arrow directions,
           pathlengths = [0.4,0.4,0.9,0.4,0.5,],
           trunklength = 1.5,
           edgecolor = '#2b8cbe',
           facecolor = '#2b8cbe' )

diagrams = sankey.finish()
for diagram in diagrams:
    for text in diagram.texts:
        text.set_fontsize(10);
plt.axis("off") 
plt.show()

Other diagram idea for the future: Sensitivity of $P_{elec}^{opt}$ to all the PTO parameters?

[jtgrasb]

@dtgaebe Thank you for following up on this! I didn't like the previous plots because they had differing input power values in the same sea state, so I think this will make a lot more sense to readers/users because the input power will be the same within the sea state.

[jtgrasb]

I've worked a bit more on this feature and wanted to provide some updates. In order to better represent a case which requires large electrical power inputs, I added an if statement to determine if the power is positive. Then, it makes the power negative, but colors the text red and manually adds a negative sign to indicate that input power is required.

The diagram is left unchanged if electrical power is harvested.

if P_elec < 0:
    Power_flows = [P_opt_excitation_6_10, P_PTO_loss, -1*P_radiated_6_4, -1*P_unused, P_elec, ]

    fig = plt.figure(figsize = [6,4])
    ax = fig.add_subplot(1, 1, 1,)
    sankey = Sankey(ax=ax, 
                    scale=.5/P_max_absorbable_6_10,
                    offset= 0,
                    format = '%.0f W',shoulder = 0.02)
    sankey.add(flows=Power_flows, 
           labels = ['$P_{e,opt}$', 
                      '$P_{loss,PTO}$', 
                      '$P_{rad}$', 
                      '$P_{unused}$', 
                      '$P_{elec}$'], 
           orientations=[0, -1, -1, -1, 0],#arrow directions,
           pathlengths = [0.2,0.1,0.3,0.1,0.2,],
           trunklength = .4,
           edgecolor = '#2b8cbe',
           facecolor = '#2b8cbe' )

    diagrams = sankey.finish()

elif P_elec > 0:
    P_elec = -abs(P_elec)
    P_tot_out = -(P_PTO_loss + -1*P_radiated_6_4 + -1*P_unused + P_elec)

    Power_flows = [P_opt_excitation_6_10, P_PTO_loss, -1*P_radiated_6_4, -1*P_unused, P_elec, ]

    fig = plt.figure(figsize = [6,4])
    ax = fig.add_subplot(1, 1, 1,)
    sankey = Sankey(ax=ax, 
                    scale=.5/P_tot_out,
                    offset= 0,
                    format = '%.0f W',shoulder = 0.02)

    sankey.add(flows=Power_flows, 
               labels = ['$P_{e,opt}$', 
                          '$P_{loss,PTO}$', 
                          '$P_{rad}$', 
                          '$P_{unused}$', 
                          '$P_{elec}$'], 
               orientations=[0, -1, -1, -1, 0],#arrow directions,
               pathlengths = [0.2,0.1,0.3,0.1,0.2,],
               trunklength = .4,
               edgecolor = '#2b8cbe',
               facecolor = '#2b8cbe' )

    diagrams = sankey.finish()

    input_text = diagrams[0].texts[-1].get_text()
    mid_text = '-'
    new_text = input_text[0:11] + mid_text + input_text[11:]
    diagrams[0].texts[-1].set_text(new_text)
    diagrams[0].texts[-1].set_color('r')

for diagram in diagrams:
    for text in diagram.texts:
        text.set_fontsize(12);

plt.axis("off") 
plt.show()

I was also looking and found this site which provides some examples of connecting Sankey diagrams to show multiple flows. I think this may be useful in the future for more detailed diagrams to better show the transition from mechanical to electrical power.

[dtgaebe]

Jantzen had a Sankey diagram up in a call with DOE and the question about the power we're leaving on the table came up. This emphasizes the need to improve the visualization that somehow communicates that the unused potential is actually undesirable, or unnecessary.

[jtgrasb]

I ran into some issues when the BEM data has negative terms for linear damping. This causes the P_opt_excitation_6_10 to be very large after check_linear_damping() adds linear friction:

In order to solve this issue, I only apply the check_linear_damping() function after calculating P_opt_excitation_6_10

P_mech = pto.mechanical_average_power(wec, x_wec, x_opt, waves, nsubsteps)
P_elec = pto.average_power(wec, x_wec, x_opt, waves, nsubsteps)

hd = wot.linear_hydrodynamics(bem_data, mass, stiffness)
#hd = wot.check_linear_damping(hd)
Zi = wot.hydrodynamic_impedance(hd)
Rad_res = np.real(Zi.squeeze())

Fex = wec_fdom.force.sel(type=['Froude_Krylov', 'diffraction']).sum('type')
Vel = wec_fdom.vel

P_max_absorbable_6_10 = (np.abs(Fex)**2/(8*Rad_res) ).squeeze().sum('omega').item()
P_opt_excitation_6_10 = 2*P_max_absorbable_6_10

hd = wot.check_linear_damping(hd)
Zi = wot.hydrodynamic_impedance(hd)
Rad_res = np.real(Zi.squeeze())

P_radiated_6_4 = ((1/2)*(Rad_res * np.abs(Vel)**2 ).squeeze().sum('omega').item())
P_excited_6_3 = (1/4)*(Fex*np.conjugate(Vel) + np.conjugate(Fex)*Vel ).squeeze().sum('omega').item()
P_PTO_absorbed_6_2 = P_excited_6_3 - P_radiated_6_4
assert(np.isclose(P_PTO_absorbed_6_2, -1*P_mech)) #assert that solver solution matches
P_PTO_loss = P_mech - P_elec
P_unused = P_opt_excitation_6_10 - P_excited_6_3

Power_flows = [P_opt_excitation_6_10, P_PTO_loss, -1*P_radiated_6_4, -1*P_unused, P_elec, ]

[dforbush2]

Jeff have you tried badBEMfix_FCM from my PR on WECSim? That should clean up the negative BEM terms

[jtgrasb]

Good point Dom, I'll try it out!