[Theory or Implementation] Extracting State-Space Matrices

[ajay-menon-iitkgp]

Hi there!

Is your question request related to a problem? Please describe.
I would like to extract specific variables from the BEMIO HDF5 output file for the RM3 device using the h5read() function. To keep things simple, I am directly using the rm3.h5 file available as an example in WEC-Sim and am only concerned with the heave mode of the float (body 1). This extraction works for most of the variables I require but yields unexpected values for the state space matrices.

Personal project context
I am extracting these variables to solve Cummins equation in Modelica using a state space approach. As matrix manipulation from a HDF5 file is rather tedious in OpenModelica, I must export the required hydro data into a '.mat' structure.

Describe the theory or implementation approach that you have a question
My understanding of WEC-Sim's implementation of the state space radiation convolution is as follows (please correct me if I am mistaken): All 4 matrices ($A$, $B$, $C$, $D$) have a maximum order of 10 for each mode. The required order is based on an R2 threshold. For example, the RM3 float (body1) achieves this threshold with an order of 2. Since heave represents the 3rd index, the state-space $A$ matrix of body 1 in the heave mode should be $A(:,:,3,3)$. After extracting this $10 \times 10$ matrix, I remove all the rows and columns that are entirely $0$ (to reduce the order to 2). Likewise, for the $B$ matrix, I extract a $10 \times 1$ column vector using $B(:,:,3,3)$ and then remove the rows that are $0$. I have implemented this using the following MATLAB code block:

filename = ""; % Insert address to rm3.h5 file

% This section works correctly - extracting hydrostatic coefficients/simulation parameters
rho = h5read(filename,'/simulation_parameters/rho');
g = h5read(filename,'/simulation_parameters/g');
Ainf33 = rho * h5read(filename,'/body1/hydro_coeffs/added_mass/inf_freq',[3 3],[1 1]);
Khs33 = rho * g * h5read(filename,'/body1/hydro_coeffs/linear_restoring_stiffness',[3 3],[1 1]);

% Read in the entire 4D dataset for each state-space matrix
ss_rad33.A = h5read(filename,'/body1/hydro_coeffs/radiation_damping/state_space/A/all');
ss_rad33.B = h5read(filename,'/body1/hydro_coeffs/radiation_damping/state_space/B/all');
ss_rad33.C = h5read(filename,'/body1/hydro_coeffs/radiation_damping/state_space/C/all');
ss_rad33.D = h5read(filename,'/body1/hydro_coeffs/radiation_damping/state_space/D/all');

% Matrix transformation to remove non-heave modes
ss_rad33.A = squeeze(hydroCoeff.ss_rad33.A(:,:,3,3));
ss_rad33.B = squeeze(hydroCoeff.ss_rad33.B(:,:,3,3));
ss_rad33.C = squeeze(hydroCoeff.ss_rad33.C(:,:,3,3));
ss_rad33.D = squeeze(hydroCoeff.ss_rad33.D(3,3));

% Strip rows and columns with all 0 values
ss_rad33.A( ~any(ss_rad33.A,2), : ) = []; ss_rad33.A( :, ~any(ss_rad33.A,1) ) = [];
ss_rad33.B( ~any(ss_rad33.B,2), : ) = []; ss_rad33.B( :, ~any(ss_rad33.B,1) ) = [];
ss_rad33.C( ~any(ss_rad33.C,2), : ) = []; ss_rad33.C( :, ~any(ss_rad33.C,1) ) = [];
ss_rad33.D( ~any(ss_rad33.D,2), : ) = []; ss_rad33.D( :, ~any(ss_rad33.D,1) ) = [];

While the coefficients are correctly extracted, the SS matrices appear to be incorrect judging from the values. For instance, comparing the extracted matrices from BEMIO with more typical ones I would expect:

$$ A_{BEMIO} = {\left\lbrack \matrix{-0.9350 & 1.0048 \cr -1.0048 & -0.0016} \right\rbrack} $$

$$ A_{typical} = {\left\lbrack \matrix{0 & 1 \cr -1.0112 & -0.9366} \right\rbrack} $$

$$ B_{BEMIO} = {\left\lbrack \matrix{-6.4077 & 0.2535} \right\rbrack} $$

$$ B_{typical} = {\left\lbrack \matrix{6.8e+05 & -5.8e+05} \right\rbrack} $$

$$ C_{BEMIO} = {\left\lbrack \matrix{-106.7948 & -4.2245} \right\rbrack} $$

$$ C_{typical} = {\left\lbrack \matrix{1 & 0} \right\rbrack} $$

$$ D_{BEMIO} = {\left\lbrack \matrix{12.8728} \right\rbrack} $$

$$ D_{typical} = {\left\lbrack \matrix{0} \right\rbrack} $$

I've obtained the typical matrices for the RM3 float using an external tool. The exported BEMIO matrices share the same dimension, but the values are all over the place. Another cause for concern is that while the $D$ matrix should be $0$ to maintain the system's causality, the value extracted from BEMIO is non-zero.

Describe the type of conclusion or resolution you are looking for
I've reviewed the WEC-Sim documentation and past issues, but nothing stands out. I used issue #162 as the basis for my understanding and code. I would appreciate your thoughts on whether the matrices extracted from BEMIO are correct (even though $D$ is non-zero) and if my MATLAB code/approach is correct.

Additional context
N/A

Cheers,
Ajay

[salhus]

Hi @ajay-menon-iitkgp,

Thanks for your observations. Some thoughts,

The matrices you are extracting are not the same as the typical matrices. For some context, refer,
https://www.sciencedirect.com/science/article/pii/0141118796000028
https://www.sciencedirect.com/science/article/pii/S0029801805000946

To your question regarding the matrix extraction, refer the lines 769-809, on how the extracted 4D matrices ought to be reshaped,

To extract the matrices, one could use the command line and exploit the protected methods in the body class shown above. Something like this, for instance,

1. Run WEC-Sim,
2. irfInfAddedMassAndDamping(body(1),50,1,1000,1),
3. body(1).hydroForce.ssRadf.A
   Let us know if reshaping the matrices as shown above resolves your issue.

Cheers,
sal

[ajay-menon-iitkgp]

@salhus, thank you for the suggestions! I will try reshaping the matrices using this method and get back to you.

To clarify: obtaining the final A, B & C matrices without running WEC-Sim requires that I implement a code similar to the image you shared that iteratively combines the intermediate matrices. Is this correct?

Cheers,
Ajay

[ajay-menon-iitkgp]

Hi @salhus. I've tried to reshape the matrices using the method shown above but am facing an issue. Consider matrix $A$ of and line no. $797$ of bodyClass.m.

Af(size(Af,1)+1:size(Af,1)+arraySize,size(Af,2)+1:size(Af,2)+arraySize) = obj.hydroData.hydro_coeffs.radiation_damping.state_space.A.all(ii,jj,1:arraySize,1:arraySize);

Variable arraySize is the optimized order of the SS matrix (stored as matrix state_space.it in the rm3.h5 file).

For specific DOFs, the order of this intermediate matrix is 10, meaning that the 4th index of A.all(ii,jj,1:arraySize,1:arraySize) in line $797$ will be $1:10$. However, the size of A.all is $10 \times 10 \times 12 \times 6$, i.e., it is bounded up to 6 elements in the 4th position. This conflict in the 4th argument generates an error in MATLAB.

I switched the arguments so that the matrices fit:

Af(...) = obj. ... .A.all(1:arraySize,1:arraySize,ii,jj);

But this yields an incorrect matrix that does not match the WEC-Sim SS output.

While I continue figuring this out, could you please help me out with 2 queries:

1. I do not understand how line $797$ works despite the conflict between 1:arraySize and the size of the 4th index of A.all.

2. This might be an elementary question: how can I generate the SS matrices for a specific DOF (say, heave) from the simulation ouptut workspace? Viewing body(1).hydroForce.ssRadf.A leaves me with a $192 \times 192$ matrix and I am unsure how to proceed.

Thank you for the assistance!

Cheers,
Ajay

[salhus]

Hi @ajay-menon-iitkgp Ajay,

Sorry I couldn't get to back to you last week. I am guessing you have BEM data and basically need a way to estimate and easily access the final SS model after the canonical forms have been processed. Here is some code I wrote a while back that you may find helpful. You could also check out the FOAMM functionality in BEM Rosetta.
compare_estimated_Falnes_Euler.zip

Let me know if this helps.

Cheers,
sal

[ajay-menon-iitkgp]

Thank you for sending the code, @salhus. I'm currently on break, so I will need a week or two to test it out. I will also take a look at the FOAMM functionality and see if that simplifies the process.

Good day!

[salhus]

Hi @ajay-menon-iitkgp ,

Hope you are doing well. This issue has been open for a bit. I was wondering if you still wanted to keep this open? If not we could close this.

Cheers,
sal

[ajay-menon-iitkgp]

Hi @salhus. Apologies for not getting back earlier. I was finally able to resolve the issue - here's a summary of the solution.

Problem: Solving Cummins equation using the SS matrices directly imported from BEMIO returned incorrect results.

Reason: The BEMIO matrices are only intermediate forms that WEC-Sim manipulates before solving Cummins equation.

Solution: The issue was solved by replicating the matrix manipulation performed by WEC-Sim's internal functions. Transpose each SS matrices and scale it by the water density $rho$. The $D$ matrix is set to $0$.

Thank you for your help with this problem. Have a great day!

Cheers,
Ajay

[salhus]

Thanks for the update @ajay-menon-iitkgp ! That makes sense. I am glad it is resolved.