Better scaling and/or initial guess

[ryancoe]

Do one or both of the following:

• Run w/o inequality constraints
• find the analytical solution (v = Fe/2/Ri)

From these, you can do one or both of the following:

• set x0
• set scaling

related to #4 #48

[ryancoe]

Note that if we use the find the analytical solution (v = Fe/2/Ri) approach as suggested above, we would ideally do this while considering the PTO impedance as being developed in #38. However, we cannot necessarily be sure how the user's PTO will look in all cases. Thus we probably want to stick with finding the analytical solution for optimal mechanical power transfer and using this for the initial guess and scaling (in the absence of input from the user). It is definitely not perfect, but probably better than nothing at all (we can/should check whether this is indeed true).

[ryancoe]

From Falnes, eq. 3.46:

[ryancoe]

@cmichelenstrofer - I see we have a function for going from a sin and cos components (i.e., our pseudo-spectral components) to a complex spectrum (real_to_complex_amplitudes), do we have anything for doing the opposite?

[cmichelenstrofer]

We do not, but it would be a good addition and easy to implement.

[ryancoe]

I've done a couple versions of this in #89. I'm going to note the successes/failures here for posterity:

• using analytic hydrodynamic solution for guess and scaling: I was only able to use this to set an initial guess for x_wec_0. Because it does not give any information on x_opt or obj, you cannot use this to set scaling and the effectiveness of setting an initial guess for only part of the solution may be limited. Nonetheless, I have added functions for the user to do this if helpful.

• run w/o inequality constraints first, then add constraints back: This appears to have mixed effectiveness. In the WaveBot case that I added to our tests, this approach does appear to improve the solution speed. However, when I tried to apply this to our examples, I did not see the same improvements.

In summary, neither of these approaches is a silver bullet, but hopefully do add some benefit to users.

[cmichelenstrofer]

@rebeccamccabe and Maha suggested using the Hessian to set the scale. They will share some code. The Hessian and scales would depend on the initial guess.

[rebeccamccabe]

My implementation of hessian scaling (in matlab) is https://github.com/symbiotic-engineering/MDOcean/blob/code-cleanup/mdocean/optimization/run_solver.m. The method came from Pappalambros & Wilde Principles of Optimal Design, section 8.3. I run the unscaled optimization for a few (8) iterations, and get the hessian from that point, use it to scale the problem, and then continue optimizing. My code gets the hessian from finite difference, but for WecOptTool you can perhaps get it from the autodiff. Scaling using the hessian as opposed to based on the magnitude of the decision variables ensures that the effect of each decision variable on the objective is comparable.

[ryancoe]

We looked into this a bit more today... a couple of notes:

Hessian

For an unstructured controller, the second derivative (diagonal of the Hessian) of our objective function (e.g., mech. power) is always zero.

$$ \frac{\delta ^2 P}{ \delta F^2} = \frac{\delta ^2 P}{ \delta v^2} \equiv 0 $$

Note, however, that mixed second order partial derivatives can be non-zero.

$$ \frac{\delta ^2 P}{ \delta v \delta F} \not\equiv 0 $$

Note also that if you were using a structured controller (e.g., PI), you would have some non-zero elements on the diagonal of Hessian.

Initial guess

While it seems possible that for specific problems you can come up with better initial guess strategies, a crude workable solution appears to be using the wave amplitude spectrum to give an initial guess. Logically, this would give you a good answer for a "wave follower" body. Initial experiments with this idea look promising.

[cmichelenstrofer]

From #225: the average mechanical power is

$$\bar{P}_{mech} = \frac{\omega _1}{2} \sum _{n=1}^{N-1}\biggl(n \left(f _{a,n} x _{b,n} - f _{b,n} x _{a,n} \right)\biggr) $$

so

$$ \frac{\partial^2 \bar{P}_{mech}}{\partial x _{a,n}^2} = 0 $$

similarly for $x _{b,n}, f _{a,n}, f _{b,n}$

and

$$\frac{\partial^2 \bar{P}_{mech}}{\partial x _{b,n}\partial f _{a,n} } = \frac{\omega _1}{2} n$$

and similar for the other cross terms.

Hessian (50 frequencies, x = [x_wec, x_opt]):

Lower diagonal:

[cmichelenstrofer]

@rebeccamccabe we found out that the diagonal of the Hessian is zero when using an unstructured controller (see examples above). What type of problem did you use the Hessian for? Was it a structured controller?

[rebeccamccabe]

My problem was not WecOptTool or pseudospectral at all, the design variables were geometric dimensions not control/dynamics harmonics, and the objective was LCOE not power, so the hessian was nonzero in most elements of the matrix.

…

On Wed, Apr 3, 2024 at 6:51 PM Carlos A. Michelén Ströfer < ***@***.***> wrote: @rebeccamccabe <https://github.com/rebeccamccabe> we found out that the diagonal of the Hessian is zero when using an unstructured controller (see examples above). What type of problem did you use the Hessian for? Was it a structured controller? — Reply to this email directly, view it on GitHub <#74 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AHUJPSXPTCDWMXEGWOILSATY3SBVZAVCNFSM5QM6W352U5DIOJSWCZC7NNSXTN2JONZXKZKDN5WW2ZLOOQ5TEMBTGU3TGOBRGUZA> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[dtgaebe]

Analog to @cmichelenstrofer's comment above, but using a PI controller the force components become a function of the dynamic states. Consequently, the Hessian does have non-zero values on it's diagonal. Let's focus on the components inside the freq array (not DC and Nyquist, $n \in [1, nfreq-1]$)

$$f_{a,n} = K_p x_{b,n} \omega_n + K_i x_{a,n}, \quad f_{b,n} = -K_p x_{a,n} \omega_n + K_{i} x_{b,n}.$$

And the mechanical power components

$$P_{a,n} = f_{a,n} \omega_n x_{b,n} , \quad P_{b,n} = f_{b,n}(- \omega x_{a,n}) $$

have non-zero values for the Hessian

$$\frac{\partial^2 P_{a,n}}{\partial x_{b,n}^2} = \frac{\partial^2 P_{b,n}}{\partial x_{a,n}^2} = K_p \omega_n^2$$

i.e. the Hessian is a quadratic function of the radial frequency scaled by the proportional gain

The numeric values for the Hessian come in pairs, when taking every other component, we can numerically confirm that the Hessian looks as we would expect.