14497012271

Committed 16 Apr 2025 03:50PM UTC coverage: 94.583%. Remained the same

Build # 14497012271

Build Type

Pull #408

github

Committed by

web-flow

Commit Message

Merge 4b947b8d0 into 1a098df54

Pull Request Pull Request #408: Sphinx build argument bug fix

Run Details

2776 of 2935 relevant lines covered (94.58%)

5.67 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

98.07

/tests/test_integration.py

""" Integration tests spanning WecOptTool.
"""
import pytest
from pytest import approx
import wecopttool as wot
import capytaine as cpy
import autograd.numpy as np
from scipy.optimize import Bounds
import xarray as xr


kplim = -1e1
min_damping = 45

@pytest.fixture(scope="module")
def f1():
    """Fundamental frequency [Hz]"""
    return 0.05


@pytest.fixture(scope="module")
def nfreq():
    """Number of frequencies in the array"""
    return 50

@pytest.fixture(scope='module')
def pto():
    """Basic PTO: unstructured, 1 DOF, mechanical power."""
    ndof = 1
    kinematics = np.eye(ndof)
    pto = wot.pto.PTO(ndof, kinematics)
    return pto


@pytest.fixture(scope='module')
def p_controller_pto():
    """Basic PTO: proportional (P) controller, 1 DOF, mechanical power."""
    ndof = 1
    pto = wot.pto.PTO(ndof=ndof, kinematics=np.eye(ndof),
                      controller=wot.pto.controller_p,
                      names=["P controller PTO"])
    return pto


@pytest.fixture(scope='module')
def pi_controller_pto():
    """Basic PTO: proportional-integral (PI) controller, 1 DOF, mechanical
    power."""
    ndof = 1
    pto = wot.pto.PTO(ndof=ndof, kinematics=np.eye(ndof),
                      controller=wot.pto.controller_pi,
                      names=["PI controller PTO"])
    return pto


@pytest.fixture(scope="module")
def fb():
    """Capytaine FloatingBody object"""
    try:
        import wecopttool.geom as geom
    except ImportError:
        pytest.skip(
            'Skipping integration tests due to missing optional geometry ' +
            'dependencies. Run `pip install wecopttool[geometry]` to run ' +
            'these tests.'
            )
    mesh_size_factor = 0.5
    wb = geom.WaveBot()
    mesh = wb.mesh(mesh_size_factor)
    fb = cpy.FloatingBody.from_meshio(mesh, name="WaveBot")
    fb.add_translation_dof(name="Heave")
    return fb


@pytest.fixture(scope="module")
def hydro_data(f1, nfreq, fb):
    """Boundary element model (Capytaine) results (with friction added)"""
    freq = wot.frequency(f1, nfreq, False)
    hydro_data = wot.run_bem(fb, freq)
    hd = wot.add_linear_friction(hydro_data)
    hd = wot.check_radiation_damping(hd, min_damping=min_damping)
    return hd


@pytest.fixture(scope='module')
def regular_wave(f1, nfreq):
    """Single frequency wave"""
    wfreq = 0.3
    wamp = 0.0625
    wphase = 0
    wdir = 0
    waves = wot.waves.regular_wave(f1, nfreq, wfreq, wamp, wphase, wdir)
    return waves


@pytest.fixture(scope='module')
def long_crested_wave(f1, nfreq):
    """Idealized (Pierson-Moskowitz) spectrum wave"""
    freq = wot.frequency(f1, nfreq, False)
    fp = 0.3
    hs = 0.0625*1.9
    spec_fun = lambda f: wot.waves.pierson_moskowitz_spectrum(freq=f, 
                                                              fp=fp, 
                                                              hs=hs)
    efth = wot.waves.omnidirectional_spectrum(f1=f1, nfreq=nfreq, 
                                              spectrum_func=spec_fun,
                                              )
    waves = wot.waves.long_crested_wave(efth, nrealizations=2)
    return waves


@pytest.fixture(scope='module')
def wec_from_bem(f1, nfreq, hydro_data, fb, pto):
    """Simple WEC: 1 DOF, no constraints."""
    f_add = {"PTO": pto.force_on_wec}
    wec = wot.WEC.from_bem(hydro_data, f_add=f_add)
    return wec


@pytest.fixture(scope='module')
def wec_from_floatingbody(f1, nfreq, fb, pto):
    """Simple WEC: 1 DOF, no constraints."""
    f_add = {"PTO": pto.force_on_wec}
    wec = wot.WEC.from_floating_body(fb, f1, nfreq, f_add=f_add, 
                                     min_damping=min_damping)
    return wec


@pytest.fixture(scope='module')
def wec_from_impedance(hydro_data, pto, fb):
    """Simple WEC: 1 DOF, no constraints."""
    bemc = hydro_data.copy()
    omega = bemc['omega'].values
    w = np.expand_dims(omega, [1,2])
    A = bemc['added_mass'].values
    B = bemc['radiation_damping'].values
    fb.center_of_mass = [0, 0, 0]
    fb.rotation_center = fb.center_of_mass
    fb = fb.copy(name=f"{fb.name}_immersed").keep_immersed_part()
    mass = bemc['inertia_matrix'].values
    hstiff = bemc['hydrostatic_stiffness'].values
    K = np.expand_dims(hstiff, 2)
    
    freqs = omega / (2 * np.pi)
    impedance = (A + mass)*(1j*w) + B + K/(1j*w)
    exc_coeff = hydro_data['Froude_Krylov_force'] + hydro_data['diffraction_force']
    f_add = {"PTO": pto.force_on_wec}

    wec = wot.WEC.from_impedance(freqs, impedance, exc_coeff, hstiff, f_add, 
                                 min_damping=min_damping)
    return wec


@pytest.fixture(scope='module')
def resonant_wave(f1, nfreq, fb, hydro_data):
    """Regular wave at natural frequency of the WEC"""
    hd = wot.add_linear_friction(hydro_data)
    Zi = wot.hydrodynamic_impedance(hd)
    wn = Zi['omega'][np.abs(Zi).argmin(dim='omega')].item()
    waves = wot.waves.regular_wave(f1, nfreq, freq=wn/2/np.pi, amplitude=0.1)
    return waves


def test_solve_callback(wec_from_bem, regular_wave, pto, nfreq, capfd):
    """Check that user can set a custom callback"""

    cbstring = 'hello world!'

    def my_callback(my_wec, x_wec, x_opt, wave):
        print(cbstring)

    _ = wec_from_bem.solve(regular_wave,
                           obj_fun=pto.average_power,
                           nstate_opt=2*nfreq,
                           scale_x_wec=1.0,
                           scale_x_opt=0.01,
                           scale_obj=1e-1,
                           callback=my_callback,
                           optim_options={'maxiter': 1})

    out, err = capfd.readouterr()

    assert out.split('\n')[0] == cbstring


@pytest.mark.parametrize("bounds_opt",
                         [Bounds(lb=kplim, ub=0), ((kplim, 0),)])
def test_solve_bounds(bounds_opt, wec_from_bem, regular_wave,
                      p_controller_pto):
    """Confirm that bounds are not violated and scale correctly when
    passing bounds argument as both as Bounds object and a tuple"""

    # replace unstructured controller with proportional controller
    wec_from_bem.forces['PTO'] = p_controller_pto.force_on_wec

    res = wec_from_bem.solve(waves=regular_wave,
                             obj_fun=p_controller_pto.average_power,
                             nstate_opt=1,
                             x_opt_0=[kplim*0.1],
                             optim_options={'maxiter': 2e1,
                                            'ftol': 1e-8},
                             bounds_opt=bounds_opt,
                             )

    assert pytest.approx(kplim, 1e-5) == res[0]['x'][-1]


def test_same_wec_init(wec_from_bem,
                       wec_from_floatingbody,
                       wec_from_impedance,
                       f1,
                       nfreq):
    """Test that different init methods for WEC class produce the same object
    """

    waves = wot.waves.regular_wave(f1, nfreq, 0.3, 0.0625)
    np.random.seed(0)
    x_wec_0 = np.random.randn(wec_from_bem.nstate_wec)
    np.random.seed(1)
    x_opt_0 = np.random.randn(wec_from_bem.nstate_wec)
    bem_res = wec_from_bem.residual(x_wec_0, x_opt_0, waves.sel(realization=0))
    fb_res = wec_from_floatingbody.residual(x_wec_0, x_opt_0, waves.sel(realization=0))
    imp_res = wec_from_impedance.residual(x_wec_0, x_opt_0, waves.sel(realization=0))
    
    assert fb_res == approx(bem_res, rel=0.01)
    assert imp_res == approx(bem_res, rel=0.01)


class TestTheoreticalPowerLimits:
    """Compare power from numerical solutions against known theoretical limits
    """

    @pytest.fixture(scope='class')
    def mass(self, fb):
        """Rigid-body mass"""
        return fb.compute_rigid_body_inertia()

    @pytest.fixture(scope='class')
    def hstiff(self, fb):
        """Hydrostatic stiffness"""
        return fb.compute_hydrostatic_stiffness()

    @pytest.fixture(scope='class')
    def hydro_impedance(self, hydro_data):
        """Intrinsic hydrodynamic impedance"""
        Zi = wot.hydrodynamic_impedance(hydro_data)
        return Zi
    
    @pytest.fixture(scope='class')
    def unstruct_wec(self,
                     hydro_data,
                     pto):
        """WaveBot WEC object with unstructured controller"""

        f_add = {"PTO": pto.force_on_wec}
        wec = wot.WEC.from_bem(hydro_data, f_add=f_add)

        return wec

    @pytest.fixture(scope='class')
    def long_crested_wave_unstruct_res(self,
                                       unstruct_wec,
                                       long_crested_wave,
                                       pto,
                                       hydro_data,
                                       nfreq):
        """Solution for an unstructured controller with multiple long crested 
        waves"""

        f_add = {"PTO": pto.force_on_wec}
        wec = wot.WEC.from_bem(hydro_data, f_add=f_add)

        res = unstruct_wec.solve(waves=long_crested_wave,
                                 obj_fun=pto.average_power,
                                 nstate_opt=2*nfreq,
                                 x_wec_0=1e-3*np.ones(wec.nstate_wec),
                                 scale_x_wec=1e1,
                                 scale_x_opt=1e-3,
                                 scale_obj=5e-2,
                                 )

        return res

    def test_p_controller_resonant_wave(self,
                                        hydro_data,
                                        resonant_wave,
                                        p_controller_pto,
                                        hydro_impedance):
        """Proportional controller should match optimum for natural resonant
        wave"""
        
        f_add = {"PTO": p_controller_pto.force_on_wec}
        wec = wot.WEC.from_bem(hydro_data, f_add=f_add)

        res = wec.solve(waves=resonant_wave,
                        obj_fun=p_controller_pto.average_power,
                        nstate_opt=1,
                        x_wec_0=1e-1*np.ones(wec.nstate_wec),
                        x_opt_0=[-1470],
                        scale_x_wec=1e2,
                        scale_x_opt=1e-3,
                        scale_obj=1e-1,
                        optim_options={'ftol': 1e-10},
                        bounds_opt=((-1*np.infty, 0),),
                        )

        power_sol = -1*res[0]['fun']

        res_fd, _ = wec.post_process(wec, res, resonant_wave, nsubsteps=1)
        Fex = res_fd[0].force.sel(
            type=['Froude_Krylov', 'diffraction']).sum('type')
        power_optimal = (np.abs(Fex)**2/8 / np.real(hydro_impedance.squeeze())
                         ).squeeze().sum('omega').item()

        assert power_sol == approx(power_optimal, rel=0.03)

    def test_pi_controller_regular_wave(self,
                                        hydro_data,
                                        regular_wave,
                                        pi_controller_pto,
                                        hydro_impedance):
        """PI controller matches optimal for any regular wave"""

        f_add = {"PTO": pi_controller_pto.force_on_wec}
        wec = wot.WEC.from_bem(hydro_data, f_add=f_add)

        res = wec.solve(waves=regular_wave,
                        obj_fun=pi_controller_pto.average_power,
                        nstate_opt=2,
                        x_wec_0=1e-1*np.ones(wec.nstate_wec),
                        x_opt_0=[-1e3, 1e4],
                        scale_x_wec=1e2,
                        scale_x_opt=1e-3,
                        scale_obj=1e-2,
                        optim_options={'maxiter': 50},
                        bounds_opt=((-1e4, 0), (0, 2e4),)
                        )

        power_sol = -1*res[0]['fun']

        res_fd, _ = wec.post_process(wec, res, regular_wave, nsubsteps=1)
        Fex = res_fd[0].force.sel(
            type=['Froude_Krylov', 'diffraction']).sum('type')
        power_optimal = (np.abs(Fex)**2/8 / np.real(hydro_impedance.squeeze())
                         ).squeeze().sum('omega').item()

        assert power_sol == approx(power_optimal, rel=1e-4)
        
    def test_unstructured_controller_long_crested_wave(self,
                                                       unstruct_wec,
                                                       long_crested_wave,
                                                       hydro_impedance,
                                                       long_crested_wave_unstruct_res,
                                                       pto):
        """Unstructured (numerical optimal) controller matches optimal for any
        irregular (long crested) wave when unconstrained"""

        power_sol = -1*long_crested_wave_unstruct_res[0]['fun']

        res_fd, _ = unstruct_wec.post_process(unstruct_wec, long_crested_wave_unstruct_res, 
                                        long_crested_wave, 
                                        nsubsteps=1)
        Fex = res_fd[0].force.sel(
            type=['Froude_Krylov', 'diffraction']).sum('type')
        power_optimal = (np.abs(Fex)**2/8 / np.real(hydro_impedance.squeeze())
                            ).squeeze().sum('omega').item()

        assert power_sol == approx(power_optimal, rel=1e-2)

    def test_unconstrained_solutions_multiple_phase_realizations(self,
                                                                 long_crested_wave_unstruct_res):
        """Solutions for average power with an unstructured controller 
        (no constraints) match for different phase realizations"""

        pow = [res['fun'] for res in long_crested_wave_unstruct_res]

        assert pow[0] == approx(pow[1], rel=1e-4)

    def test_saturated_pi_controller(self,
                                    hydro_data,
                                    regular_wave,
                                    pto,
                                    nfreq):
        """Saturated PI controller matches constrained unstructured controller
        for a regular wave
        """

        pto_tmp = pto
        pto = {}
        wec = {}
        nstate_opt = {}
        
        # Constraint
        f_max = 2000.0

        nstate_opt['us'] = 2*nfreq
        pto['us'] = pto_tmp
        def const_f_pto(wec, x_wec, x_opt, waves):
            f = pto['us'].force_on_wec(wec, x_wec, x_opt, waves, 
                                       nsubsteps=4)
            return f_max - np.abs(f.flatten())
        wec['us'] = wot.WEC.from_bem(hydro_data,
                                     f_add={"PTO": pto['us'].force_on_wec},
                                     constraints=[{'type': 'ineq',
                                                   'fun': const_f_pto, }])
        
        
        ndof = 1
        nstate_opt['pi'] = 2
        def saturated_pi(pto, wec, x_wec, x_opt, waves=None, nsubsteps=1):
            return wot.pto.controller_pi(pto, wec, x_wec, x_opt, waves, 
                                         nsubsteps, 
                                         saturation=[-f_max, f_max])
        pto['pi'] = wot.pto.PTO(ndof=ndof,
                                kinematics=np.eye(ndof),
                                controller=saturated_pi,)
        wec['pi'] = wot.WEC.from_bem(hydro_data,
                                     f_add={"PTO": pto['pi'].force_on_wec},
                                     constraints=[])
        
        x_opt_0 = {'us': np.ones(nstate_opt['us'])*0.1,
                   'pi': [-1e3, 1e4]}
        scale_x_wec = {'us': 1e1,
                       'pi': 1e1}
        scale_x_opt = {'us': 1e-3,
                       'pi': 1e-3}
        scale_obj = {'us': 1e-2,
                     'pi': 1e-2}
        bounds_opt = {'us': None,
                      'pi': ((-1e4, 0), (0, 2e4),)}
        
        res = {}
        pto_fdom = {}
        pto_tdom = {}
        for key in wec.keys():
            res[key] = wec[key].solve(waves=regular_wave,
                            obj_fun=pto[key].average_power,
                            nstate_opt=nstate_opt[key],
                            x_wec_0=1e-1*np.ones(wec[key].nstate_wec),
                            x_opt_0=x_opt_0[key],
                            scale_x_wec=scale_x_wec[key],
                            scale_x_opt=scale_x_opt[key],
                            scale_obj=scale_obj[key],
                            optim_options={'maxiter': 200},
                            bounds_opt=bounds_opt[key]
                            )
            
            nsubstep_postprocess = 4
            pto_fdom[key], pto_tdom[key] = pto[key].post_process(wec[key], 
                                                                 res[key], 
                                                                 regular_wave, 
                                                                 nsubstep_postprocess)
        
        xr.testing.assert_allclose(pto_tdom['pi'][0].power.squeeze().mean('time'), 
                                   pto_tdom['us'][0].power.squeeze().mean('time'),
                                   rtol=1e-1)
        
        xr.testing.assert_allclose(pto_tdom['us'][0].force.max(),
                                   pto_tdom['pi'][0].force.max())

1	""" Integration tests spanning WecOptTool.
2	"""
3	import pytest	6✔
4	from pytest import approx	6✔
5	import wecopttool as wot	6✔
6	import capytaine as cpy	6✔
7	import autograd.numpy as np	6✔
8	from scipy.optimize import Bounds	6✔
9	import xarray as xr	6✔
10
11
12	kplim = -1e1	6✔
13	min_damping = 45	6✔
14
15	@pytest.fixture(scope="module")	6✔
16	def f1():	6✔
17	"""Fundamental frequency [Hz]"""
18	return 0.05	6✔
19
20
21	@pytest.fixture(scope="module")	6✔
22	def nfreq():	6✔
23	"""Number of frequencies in the array"""
24	return 50	6✔
25
26	@pytest.fixture(scope='module')	6✔
27	def pto():	6✔
28	"""Basic PTO: unstructured, 1 DOF, mechanical power."""
29	ndof = 1	6✔
30	kinematics = np.eye(ndof)	6✔
31	pto = wot.pto.PTO(ndof, kinematics)	6✔
32	return pto	6✔
33
34
35	@pytest.fixture(scope='module')	6✔
36	def p_controller_pto():	6✔
37	"""Basic PTO: proportional (P) controller, 1 DOF, mechanical power."""
38	ndof = 1	6✔
39	pto = wot.pto.PTO(ndof=ndof, kinematics=np.eye(ndof),	6✔
40	controller=wot.pto.controller_p,
41	names=["P controller PTO"])
42	return pto	6✔
43
44
45	@pytest.fixture(scope='module')	6✔
46	def pi_controller_pto():	6✔
47	"""Basic PTO: proportional-integral (PI) controller, 1 DOF, mechanical
48	power."""
49	ndof = 1	6✔
50	pto = wot.pto.PTO(ndof=ndof, kinematics=np.eye(ndof),	6✔
51	controller=wot.pto.controller_pi,
52	names=["PI controller PTO"])
53	return pto	6✔
54
55
56	@pytest.fixture(scope="module")	6✔
57	def fb():	6✔
58	"""Capytaine FloatingBody object"""
59	try:	6✔
60	import wecopttool.geom as geom	6✔
61	except ImportError:	×
62	pytest.skip(	×
63	'Skipping integration tests due to missing optional geometry ' +
64	'dependencies. Run `pip install wecopttool[geometry]` to run ' +
65	'these tests.'
66	)
67	mesh_size_factor = 0.5	6✔
68	wb = geom.WaveBot()	6✔
69	mesh = wb.mesh(mesh_size_factor)	6✔
70	fb = cpy.FloatingBody.from_meshio(mesh, name="WaveBot")	6✔
71	fb.add_translation_dof(name="Heave")	6✔
72	return fb	6✔
73
74
75	@pytest.fixture(scope="module")	6✔
76	def hydro_data(f1, nfreq, fb):	6✔
77	"""Boundary element model (Capytaine) results (with friction added)"""
78	freq = wot.frequency(f1, nfreq, False)	6✔
79	hydro_data = wot.run_bem(fb, freq)	6✔
80	hd = wot.add_linear_friction(hydro_data)	6✔
81	hd = wot.check_radiation_damping(hd, min_damping=min_damping)	6✔
82	return hd	6✔
83
84
85	@pytest.fixture(scope='module')	6✔
86	def regular_wave(f1, nfreq):	6✔
87	"""Single frequency wave"""
88	wfreq = 0.3	6✔
89	wamp = 0.0625	6✔
90	wphase = 0	6✔
91	wdir = 0	6✔
92	waves = wot.waves.regular_wave(f1, nfreq, wfreq, wamp, wphase, wdir)	6✔
93	return waves	6✔
94
95
96	@pytest.fixture(scope='module')	6✔
97	def long_crested_wave(f1, nfreq):	6✔
98	"""Idealized (Pierson-Moskowitz) spectrum wave"""
99	freq = wot.frequency(f1, nfreq, False)	6✔
100	fp = 0.3	6✔
101	hs = 0.0625*1.9	6✔
102	spec_fun = lambda f: wot.waves.pierson_moskowitz_spectrum(freq=f,	6✔
103	fp=fp,
104	hs=hs)
105	efth = wot.waves.omnidirectional_spectrum(f1=f1, nfreq=nfreq,	6✔
106	spectrum_func=spec_fun,
107	)
108	waves = wot.waves.long_crested_wave(efth, nrealizations=2)	6✔
109	return waves	6✔
110
111
112	@pytest.fixture(scope='module')	6✔
113	def wec_from_bem(f1, nfreq, hydro_data, fb, pto):	6✔
114	"""Simple WEC: 1 DOF, no constraints."""
115	f_add = {"PTO": pto.force_on_wec}	6✔
116	wec = wot.WEC.from_bem(hydro_data, f_add=f_add)	6✔
117	return wec	6✔
118
119
120	@pytest.fixture(scope='module')	6✔
121	def wec_from_floatingbody(f1, nfreq, fb, pto):	6✔
122	"""Simple WEC: 1 DOF, no constraints."""
123	f_add = {"PTO": pto.force_on_wec}	6✔
124	wec = wot.WEC.from_floating_body(fb, f1, nfreq, f_add=f_add,	6✔
125	min_damping=min_damping)
126	return wec	6✔
127
128
129	@pytest.fixture(scope='module')	6✔
130	def wec_from_impedance(hydro_data, pto, fb):	6✔
131	"""Simple WEC: 1 DOF, no constraints."""
132	bemc = hydro_data.copy()	6✔
133	omega = bemc['omega'].values	6✔
134	w = np.expand_dims(omega, [1,2])	6✔
135	A = bemc['added_mass'].values	6✔
136	B = bemc['radiation_damping'].values	6✔
137	fb.center_of_mass = [0, 0, 0]	6✔
138	fb.rotation_center = fb.center_of_mass	6✔
139	fb = fb.copy(name=f"{fb.name}_immersed").keep_immersed_part()	6✔
140	mass = bemc['inertia_matrix'].values	6✔
141	hstiff = bemc['hydrostatic_stiffness'].values	6✔
142	K = np.expand_dims(hstiff, 2)	6✔
143
144	freqs = omega / (2 * np.pi)	6✔
145	impedance = (A + mass)(1jw) + B + K/(1j*w)	6✔
146	exc_coeff = hydro_data['Froude_Krylov_force'] + hydro_data['diffraction_force']	6✔
147	f_add = {"PTO": pto.force_on_wec}	6✔
148
149	wec = wot.WEC.from_impedance(freqs, impedance, exc_coeff, hstiff, f_add,	6✔
150	min_damping=min_damping)
151	return wec	6✔
152
153
154	@pytest.fixture(scope='module')	6✔
155	def resonant_wave(f1, nfreq, fb, hydro_data):	6✔
156	"""Regular wave at natural frequency of the WEC"""
157	hd = wot.add_linear_friction(hydro_data)	6✔
158	Zi = wot.hydrodynamic_impedance(hd)	6✔
159	wn = Zi['omega'][np.abs(Zi).argmin(dim='omega')].item()	6✔
160	waves = wot.waves.regular_wave(f1, nfreq, freq=wn/2/np.pi, amplitude=0.1)	6✔
161	return waves	6✔
162
163
164	def test_solve_callback(wec_from_bem, regular_wave, pto, nfreq, capfd):	6✔
165	"""Check that user can set a custom callback"""
166
167	cbstring = 'hello world!'	6✔
168
169	def my_callback(my_wec, x_wec, x_opt, wave):	6✔
170	print(cbstring)	6✔
171
172	_ = wec_from_bem.solve(regular_wave,	6✔
173	obj_fun=pto.average_power,
174	nstate_opt=2*nfreq,
175	scale_x_wec=1.0,
176	scale_x_opt=0.01,
177	scale_obj=1e-1,
178	callback=my_callback,
179	optim_options={'maxiter': 1})
180
181	out, err = capfd.readouterr()	6✔
182
183	assert out.split('\n')[0] == cbstring	6✔
184
185
186	@pytest.mark.parametrize("bounds_opt",	6✔
187	[Bounds(lb=kplim, ub=0), ((kplim, 0),)])
188	def test_solve_bounds(bounds_opt, wec_from_bem, regular_wave,	6✔
189	p_controller_pto):
190	"""Confirm that bounds are not violated and scale correctly when
191	passing bounds argument as both as Bounds object and a tuple"""
192
193	# replace unstructured controller with proportional controller
194	wec_from_bem.forces['PTO'] = p_controller_pto.force_on_wec	6✔
195
196	res = wec_from_bem.solve(waves=regular_wave,	6✔
197	obj_fun=p_controller_pto.average_power,
198	nstate_opt=1,
199	x_opt_0=[kplim*0.1],
200	optim_options={'maxiter': 2e1,
201	'ftol': 1e-8},
202	bounds_opt=bounds_opt,
203	)
204
205	assert pytest.approx(kplim, 1e-5) == res[0]['x'][-1]	6✔
206
207
208	def test_same_wec_init(wec_from_bem,	6✔
209	wec_from_floatingbody,
210	wec_from_impedance,
211	f1,
212	nfreq):
213	"""Test that different init methods for WEC class produce the same object
214	"""
215
216	waves = wot.waves.regular_wave(f1, nfreq, 0.3, 0.0625)	6✔
217	np.random.seed(0)	6✔
218	x_wec_0 = np.random.randn(wec_from_bem.nstate_wec)	6✔
219	np.random.seed(1)	6✔
220	x_opt_0 = np.random.randn(wec_from_bem.nstate_wec)	6✔
221	bem_res = wec_from_bem.residual(x_wec_0, x_opt_0, waves.sel(realization=0))	6✔
222	fb_res = wec_from_floatingbody.residual(x_wec_0, x_opt_0, waves.sel(realization=0))	6✔
223	imp_res = wec_from_impedance.residual(x_wec_0, x_opt_0, waves.sel(realization=0))	6✔
224
225	assert fb_res == approx(bem_res, rel=0.01)	6✔
226	assert imp_res == approx(bem_res, rel=0.01)	6✔
227
228
229	class TestTheoreticalPowerLimits:	6✔
230	"""Compare power from numerical solutions against known theoretical limits
231	"""
232
233	@pytest.fixture(scope='class')	6✔
234	def mass(self, fb):	6✔
235	"""Rigid-body mass"""
236	return fb.compute_rigid_body_inertia()	×
237
238	@pytest.fixture(scope='class')	6✔
239	def hstiff(self, fb):	6✔
240	"""Hydrostatic stiffness"""
241	return fb.compute_hydrostatic_stiffness()	×
242
243	@pytest.fixture(scope='class')	6✔
244	def hydro_impedance(self, hydro_data):	6✔
245	"""Intrinsic hydrodynamic impedance"""
246	Zi = wot.hydrodynamic_impedance(hydro_data)	6✔
247	return Zi	6✔
248
249	@pytest.fixture(scope='class')	6✔
250	def unstruct_wec(self,	6✔
251	hydro_data,
252	pto):
253	"""WaveBot WEC object with unstructured controller"""
254
255	f_add = {"PTO": pto.force_on_wec}	6✔
256	wec = wot.WEC.from_bem(hydro_data, f_add=f_add)	6✔
257
258	return wec	6✔
259
260	@pytest.fixture(scope='class')	6✔
261	def long_crested_wave_unstruct_res(self,	6✔
262	unstruct_wec,
263	long_crested_wave,
264	pto,
265	hydro_data,
266	nfreq):
267	"""Solution for an unstructured controller with multiple long crested
268	waves"""
269
270	f_add = {"PTO": pto.force_on_wec}	6✔
271	wec = wot.WEC.from_bem(hydro_data, f_add=f_add)	6✔
272
273	res = unstruct_wec.solve(waves=long_crested_wave,	6✔
274	obj_fun=pto.average_power,
275	nstate_opt=2*nfreq,
276	x_wec_0=1e-3*np.ones(wec.nstate_wec),
277	scale_x_wec=1e1,
278	scale_x_opt=1e-3,
279	scale_obj=5e-2,
280	)
281
282	return res	6✔
283
284	def test_p_controller_resonant_wave(self,	6✔
285	hydro_data,
286	resonant_wave,
287	p_controller_pto,
288	hydro_impedance):
289	"""Proportional controller should match optimum for natural resonant
290	wave"""
291
292	f_add = {"PTO": p_controller_pto.force_on_wec}	6✔
293	wec = wot.WEC.from_bem(hydro_data, f_add=f_add)	6✔
294
295	res = wec.solve(waves=resonant_wave,	6✔
296	obj_fun=p_controller_pto.average_power,
297	nstate_opt=1,
298	x_wec_0=1e-1*np.ones(wec.nstate_wec),
299	x_opt_0=[-1470],
300	scale_x_wec=1e2,
301	scale_x_opt=1e-3,
302	scale_obj=1e-1,
303	optim_options={'ftol': 1e-10},
304	bounds_opt=((-1*np.infty, 0),),
305	)
306
307	power_sol = -1*res[0]['fun']	6✔
308
309	res_fd, _ = wec.post_process(wec, res, resonant_wave, nsubsteps=1)	6✔
310	Fex = res_fd[0].force.sel(	6✔
311	type=['Froude_Krylov', 'diffraction']).sum('type')
312	power_optimal = (np.abs(Fex)**2/8 / np.real(hydro_impedance.squeeze())	6✔
313	).squeeze().sum('omega').item()
314
315	assert power_sol == approx(power_optimal, rel=0.03)	6✔
316
317	def test_pi_controller_regular_wave(self,	6✔
318	hydro_data,
319	regular_wave,
320	pi_controller_pto,
321	hydro_impedance):
322	"""PI controller matches optimal for any regular wave"""
323
324	f_add = {"PTO": pi_controller_pto.force_on_wec}	6✔
325	wec = wot.WEC.from_bem(hydro_data, f_add=f_add)	6✔
326
327	res = wec.solve(waves=regular_wave,	6✔
328	obj_fun=pi_controller_pto.average_power,
329	nstate_opt=2,
330	x_wec_0=1e-1*np.ones(wec.nstate_wec),
331	x_opt_0=[-1e3, 1e4],
332	scale_x_wec=1e2,
333	scale_x_opt=1e-3,
334	scale_obj=1e-2,
335	optim_options={'maxiter': 50},
336	bounds_opt=((-1e4, 0), (0, 2e4),)
337	)
338
339	power_sol = -1*res[0]['fun']	6✔
340
341	res_fd, _ = wec.post_process(wec, res, regular_wave, nsubsteps=1)	6✔
342	Fex = res_fd[0].force.sel(	6✔
343	type=['Froude_Krylov', 'diffraction']).sum('type')
344	power_optimal = (np.abs(Fex)**2/8 / np.real(hydro_impedance.squeeze())	6✔
345	).squeeze().sum('omega').item()
346
347	assert power_sol == approx(power_optimal, rel=1e-4)	6✔
348
349	def test_unstructured_controller_long_crested_wave(self,	6✔
350	unstruct_wec,
351	long_crested_wave,
352	hydro_impedance,
353	long_crested_wave_unstruct_res,
354	pto):
355	"""Unstructured (numerical optimal) controller matches optimal for any
356	irregular (long crested) wave when unconstrained"""
357
358	power_sol = -1*long_crested_wave_unstruct_res[0]['fun']	6✔
359
360	res_fd, _ = unstruct_wec.post_process(unstruct_wec, long_crested_wave_unstruct_res,	6✔
361	long_crested_wave,
362	nsubsteps=1)
363	Fex = res_fd[0].force.sel(	6✔
364	type=['Froude_Krylov', 'diffraction']).sum('type')
365	power_optimal = (np.abs(Fex)**2/8 / np.real(hydro_impedance.squeeze())	6✔
366	).squeeze().sum('omega').item()
367
368	assert power_sol == approx(power_optimal, rel=1e-2)	6✔
369
370	def test_unconstrained_solutions_multiple_phase_realizations(self,	6✔
371	long_crested_wave_unstruct_res):
372	"""Solutions for average power with an unstructured controller
373	(no constraints) match for different phase realizations"""
374
375	pow = [res['fun'] for res in long_crested_wave_unstruct_res]	6✔
376
377	assert pow[0] == approx(pow[1], rel=1e-4)	6✔
378
379	def test_saturated_pi_controller(self,	6✔
380	hydro_data,
381	regular_wave,
382	pto,
383	nfreq):
384	"""Saturated PI controller matches constrained unstructured controller
385	for a regular wave
386	"""
387
388	pto_tmp = pto	6✔
389	pto = {}	6✔
390	wec = {}	6✔
391	nstate_opt = {}	6✔
392
393	# Constraint
394	f_max = 2000.0	6✔
395
396	nstate_opt['us'] = 2*nfreq	6✔
397	pto['us'] = pto_tmp	6✔
398	def const_f_pto(wec, x_wec, x_opt, waves):	6✔
399	f = pto['us'].force_on_wec(wec, x_wec, x_opt, waves,	6✔
400	nsubsteps=4)
401	return f_max - np.abs(f.flatten())	6✔
402	wec['us'] = wot.WEC.from_bem(hydro_data,	6✔
403	f_add={"PTO": pto['us'].force_on_wec},
404	constraints=[{'type': 'ineq',
405	'fun': const_f_pto, }])
406
407
408	ndof = 1	6✔
409	nstate_opt['pi'] = 2	6✔
410	def saturated_pi(pto, wec, x_wec, x_opt, waves=None, nsubsteps=1):	6✔
411	return wot.pto.controller_pi(pto, wec, x_wec, x_opt, waves,	6✔
412	nsubsteps,
413	saturation=[-f_max, f_max])
414	pto['pi'] = wot.pto.PTO(ndof=ndof,	6✔
415	kinematics=np.eye(ndof),
416	controller=saturated_pi,)
417	wec['pi'] = wot.WEC.from_bem(hydro_data,	6✔
418	f_add={"PTO": pto['pi'].force_on_wec},
419	constraints=[])
420
421	x_opt_0 = {'us': np.ones(nstate_opt['us'])*0.1,	6✔
422	'pi': [-1e3, 1e4]}
423	scale_x_wec = {'us': 1e1,	6✔
424	'pi': 1e1}
425	scale_x_opt = {'us': 1e-3,	6✔
426	'pi': 1e-3}
427	scale_obj = {'us': 1e-2,	6✔
428	'pi': 1e-2}
429	bounds_opt = {'us': None,	6✔
430	'pi': ((-1e4, 0), (0, 2e4),)}
431
432	res = {}	6✔
433	pto_fdom = {}	6✔
434	pto_tdom = {}	6✔
435	for key in wec.keys():	6✔
436	res[key] = wec[key].solve(waves=regular_wave,	6✔
437	obj_fun=pto[key].average_power,
438	nstate_opt=nstate_opt[key],
439	x_wec_0=1e-1*np.ones(wec[key].nstate_wec),
440	x_opt_0=x_opt_0[key],
441	scale_x_wec=scale_x_wec[key],
442	scale_x_opt=scale_x_opt[key],
443	scale_obj=scale_obj[key],
444	optim_options={'maxiter': 200},
445	bounds_opt=bounds_opt[key]
446	)
447
448	nsubstep_postprocess = 4	6✔
449	pto_fdom[key], pto_tdom[key] = pto[key].post_process(wec[key],	6✔
450	res[key],
451	regular_wave,
452	nsubstep_postprocess)
453
454	xr.testing.assert_allclose(pto_tdom['pi'][0].power.squeeze().mean('time'),	6✔
455	pto_tdom['us'][0].power.squeeze().mean('time'),
456	rtol=1e-1)
457
458	xr.testing.assert_allclose(pto_tdom['us'][0].force.max(),	6✔
459	pto_tdom['pi'][0].force.max())

sandialabs / WecOptTool / 14497012271

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous