Senthuran Ravinthrakumar
email: senthuran.ravinthrakumar@ntnu.no

Trygve Kristiansen

Bernard Molin

Babak Ommani

Coupled vessel and moonpool responses in regular and irregular waves

Keywords:

Hydrodynamic coupling between the ship and moonpool responses is investigated. Dedicated experiments are carried out in the Ocean Basin at Sintef Ocean. Three different moonpool sizes are studied, where the moonpool length is 1/20, 1/10 and 1/2 of the ship length. The two former are square moonpools, while the latter is a rectangular moonpool with a width 1/2 of the ship's beam. Numerical simulations are carried out using WAMIT. The coupling between the ship and moonpool is seen to be highly dependent on the volume of the moonpool relative to the submerged ship volume. WAMIT greatly over-predicts the moonpool responses in the proximity of resonance, which indicates that wave radiation damping is small, and that damping due to flow separation at sharp edges at the moonpool entrance is dominant. Two main nonlinear effects are observed; Swirling-type sloshing and secondary resonance. These nonlinear effects are excited in a range of wave periods where it is expected that a realistic sea environment will contain significant wave energy. This is most prominent in the largest moonpool. The piston mode shape in freely floating conditions is signficantly different relative to the one in forced heave.

Introduction

A moonpool is a vertical opening in a ship structure, which introduces an alternative option to lower objects into the sea. Another possible alternative is to off-load the object on the side of the ship. On the contrary to the alternative options, the moonpool shelters the objects from environmental loads to some extent. The moonpool, however, might exhibit significant resonant motions, which serve as a limiting factor to perform a certain operation through the moonpool. Traditionally, moonpool resonance is studied in forced heave conditions. In the present study, we will focus on resonant motions in the moonpool, and its effect on the vessel motions, in freely floating conditions. Both the piston and sloshing modes are covered. The piston mode period is not well separated from the first longitudinal sloshing mode for the largest moonpool.

When resonant motion in the moonpool is excited, the oscillations can reach several multiple of the incident wave height. In such circumstances, it might be challenging to work and perform an operation. Moreover, if the moonpool motions becomes vigorous it can lead to damages on equipment and water on the working deck. Resonant moonpool motions can be associated with large waiting time for the moonpool to calm down, which affects the operability of a vessel.

The large-amplitude moonpool oscillations in the proximity of the piston and sloshing modes are of concern with respect to large hydrodynamic loads, including slamming loads, on the structure being lowered through the moonpool. Sloshing can cause relatively large horizontal and vertical motions in the moonpool. Large horizontal moonpool oscillations can be of concern if the structure being lowered through the moonpool is vulnerable to horizontal loads. [START_REF] Molin | On the piston and sloshing modes in moonpools[END_REF] developed an analytical method to predict the natural periods and corresponding modes of the resonant motions in a moonpool by using a domain-decomposition method, assuming linear potential flow theory. The methods involve solving an eigenvalue problem, and do not predict the amplitude of the moonpool response.

Experimental and numerical investigations of the moonpool response were carried out by [START_REF] Faltinsen | Two-dimensional resonant pistonlike sloshing in a moonpool[END_REF][START_REF] Kristiansen | Application of a vortex tracking method to the piston-like behaviour in a semi-entrained vertical gap[END_REF][START_REF] Kristiansen | Gap resonance analyzed by a new domain-decomposition method combining potential and viscous flow[END_REF], where they found that potential flow theory significantly over-predicted the moonpool responses at resonance. The studies show that there are two main sources of damping of the moonpool oscillations; Radiation and viscous damping. The former is associated with radiated waves due to the communication with the outer domain, and is in general small. The latter is associated with damping due to flow separation at the moonpool entrance, and is seen to be dominant in the proximity of the piston mode. [START_REF] Ravinthrakumar | A two-dimensional numerical and experimental study of piston and sloshing resonance in moonpools with recess[END_REF] studied moonpools with recess in a two-dimensional setting, where they found dominant damping due to flow separation also in the proximity of the first sloshing mode. They discussed that the radiation damping was significantly smaller in the proximity of the first sloshing mode than around the piston mode.

In freely floating conditions, the resonant moonpool responses are different from those in forced motion tests. The resonant piston mode responses are typically considerably smaller. Also, the piston mode resonance period is changed. The latter is due to that poles in the radiation and diffraction potentials cancel out each other [START_REF] Mciver | Complex resonances in the water-wave problem for a floating structure[END_REF] at the piston mode period as predicted in e.g. forced heave. Further, if the moonpool size is large relative to the vessel, the vessel responses are influenced by mutual interaction. This was also investigated by [START_REF] Fredriksen | Wave-induced response of a floating two-dimensional body with a moonpool[END_REF]. They studied the importance of the hydrodynamic interaction effect between the moonpool and vessel motions in a two-dimensional setting. They discussed the significant coupling effects between the heave and piston mode motions. They observed that the moonpool influenced the heave responses in the vicinity of the piston mode, and that the piston mode period was shifted due to the coupling between heave and moonpool motions. The moonpool volume, V MP , was significant relative to the submerged ship volume, V S , with = V V / 0.25

MP S
. [START_REF] Ravinthrakumar | On the hydrodynamic interaction between ship and free-surface motions on vessels with moonpools[END_REF] studied the importance of the piston mode response on the heave and pitch motions in a two-dimensional and quasi three-dimensional setting with = V V / 0.022, 0.045, 0.094 and 0.20,

MP S

where they confirmed that the hydrodynamic coupling effect was clearly dependent on the V MP /V S -ratio. [START_REF] Guo | Resonant water motions within a recessing type moonpool in a drilling vessel[END_REF] performed experimental studies in an offshore basin, and investigated the resonant moonpool responses in large moonpools with recess in regular waves. They observed higher harmonic moonpool responses when sloshing occurred. We clearly distinguish between higher harmonic responses and higher modes. The former refers to the nth harmonic responses oscillating with frequencies of ω, 2ω, , nω, where ω is the circular wave frequency. The latter refers to resonant modes that are activated due to secondary resonance [2]. [START_REF] Van't Veer | Added resistance of moonpools in calm water[END_REF] presented an overview of several drillships with their main dimensions. Based on their overview, it can be found that the typical V MP / V S -ratio is around 1-6%. Over the recent years, the moonpool dimensions have increased in size, and there is a need for understanding the dominant physical aspects when the moonpool dimensions become large. This paper addresses resonant motion of the moonpool, and its effect on the ship responses. We do this by studying a ship with three different moonpool sizes; = V V / 0.015, 0.064 and 0.43

MP S
. The smallest moonpool has the dimensions of a moonpool used to typically lower Remotely Operated Vehicles (ROVs) and for recovery of divers. The second moonpool we study is moderate, and is more commonly found on drillships. The small and moderate sized moonpools are located in the fore part of the vessel, since this is common when there is a desire to utilize the deck space for other purposes. These models are inspired by an existing ship which was investigated by Sintef Ocean in 2018. The largest moonpool is significantly larger than the two former, and a moonpool with these dimensions is not built to the authors' knowledge. The purpose of studying the large moonpool is to investigate the resonant motions that can be induced in a real sea state, and the consequences these resonant modes will have on the vessel motions. [START_REF] Maisondieu | Flow dynamics in a moon-pool experimental and numerical assessment[END_REF] studied a Wellhead barge floating unit, where the moonpool length was approximately 45% of the overall length. The concept was dedicated to oil productions in deep waters. They observed significant coupling between the ship and moonpool responses. The length of the largest moonpool in the present study is 50% of the overall length.

In the present study, we address the piston and sloshing modes in a relevant range of wave periods for a ship in operational conditions in the North Sea, and discuss how these influence the vessel responses. A ship model with three different moonpool sizes is investigated by performing model tests in the Ocean Basin at Sintef Ocean. Main findings in the present study include the importance of the V MP /V S -ratio on the hydrodynamic interaction effect between the moonpool and ship responses. The importance of considering the sloshing modes for moderate and large moonpools is also addressed, especially due to excitation of higher modes due to nonlinear effects. The piston mode shape for the largest moonpool configuration in freely floating conditions is also discussed. The importance of damping due to flow separation at the moonpool entrance for the piston and sloshing modes is discussed.

The paper is organized as follows: The experimental set-up is discussed in Section 2. In Section 3 the numerical simulations are presented, where we also present a convergence study. Finally, Section 4 presents and discusses the results from the experimental and numerical studies, followed by conclusions in Section 5.

Experimental set-up

Dedicated experimental studies were carried out in the Ocean Basin at Sintef Ocean in Trondheim, Norway. The basin is 80 m long, 50 m wide and the water depth was 5 m. The basin is equipped with a hinged double-flap wavemaker. In the present study, vessel and moonpool responses are studied in long-crested waves with a ship model in operational conditions.

The model length was 4 m, and a scale of 1:34.5 is imagined. Other main particulars of the ship model are summarized in Table 1. The model was built for the present study, where the model was made such that different moonpool sizes could be tested. Photos of the model are presented in Fig. 1. The hull has verticall sides and horizontal bottom. The model is not equipped with bilge keels. We believe there are no appreciable scaling effects between model scale and full scale, since the separation point is fixed and the same in both model and full scale, due to that the present models have sharp corners.

Three different moonpool sizes are tested, referred to as MP1, MP2 and MP3 (cf. Table 1). A cradle test was conducted to determine the radius of gyration in pitch and center of mass of the naked and fully equipped models with added weights. The total mass was also measured during this process. The results were compared to calculations using the CAD software Rhino, and differences less than 0.5% were obtained. Weights in the model were placed such that zero trim and heel angles were obtained in calm water. In the present study, the origin of the coordinate system is placed in the center of the model at the mean waterline (i.e. L/2 from the aft of the model). All results in this paper are presented in model scale, unless otherwise stated.

The model was equipped with body-fixed wave probes (WPs) evenly and pairwise distributed in the moonpool, with four, six and twelve wave probes for MP1, MP2 and MP3, respectively (cf. Fig. 2). In addition, six Earth-fixed wave probes were used outside the model. Oqus position system was used to measure the motion in six directions. The model was equipped with three body-fixed accelerometers, each measuring accelerations in x-, y-and z-directions. Moreover, the model was equipped with four force rings, one in each mooring line, and two force transducers (one on each end of the model. The force measurements are not presented herein.

The vessel was moored with four horizontal mooring lines, each with a stiffness of 90 N/m. The pre-tension in the mooring line was prescribed to be 45 N. The pre-tension was checked during the model tests with the force transducers, and the measurements show that the prescribed pre-tension was achieved with a difference less than 3%.

The experiments were carried out with the model freely floating in incident waves. The tests were conducted in both regular and irregular waves. The duration of the former was three minutes, with a waiting time of six minutes between each test for the basin to calm down. The irregular waves consisted of time-series of 15 minutes, with a waiting time of 15 min between each test. The first three minutes of the irregular wave tests are disregarded during post-processing, such that we analyze a time-window where all of the wave components have reached the model.

Waves

Regular waves in the present study are generated given wave steepness and wave period. A regular wave with wave steepness, H/λ, and wave period, T, is

= x t H t kx ( , ) 2 sin( ) (1)
where k is the wave number, = k 2 / is the wavelength, = T 2 / is the circular frequency and H is the wave height.

We generate irregular waves in an attempt to be able to directly compare the Response Amplitude Operator (RAO) from the irregular wave tests with RAOs from the regular wave tests. We established a formulation for the wave spectrum, S ζ , based on the spectrum value for a regular wave given wave steepness and wave period. The irregular waves tested in the Ocean Basin are summarized in Table 2, referred to as Pink Noise for Water Waves spectrum in this paper. The specified and measured wave spectra are compared in Fig. 3. We found satisfactory agreement between the desired and measured wave spectra in the Ocean Basin, although achieved is up to 15% higher than specified for periods close to 2 s (ω ≈ 3 rad/s).

The RAO of a variable R is defined in the usual sense as

= H S S ( ) ( ) ( ) R R (2)
where H R and S R are the RAO and spectrum of the variable R, respectively. Here, we assume that the spectrum represents a stationary process and that there is a linear relationship between the wave and response spectra.

Error discussions

The model tests were conducted carefully to minimize inaccuracies and errors in measurements. We have assessed possible bias and random errors to discuss the quality of the experimental results.

The measured and specified wave heights in regular waves were compared, where satisfactory agreement between specified and measured wave heights is observed. The maximum difference in the desired and measured wave amplitude is 7%. Throughout the rest of the paper, we refer to the measured wave height. The duration of each regular wave test was three minutes, with a linear ramp up of ten seconds at the start of each test. Wave documentation tests were performed without the model, at the mean model position.

The rigid body motions as measured by Oqus position measurements were compared to acceleration measurements by integrating twice in time using standard Fourier transformation. The results showed satisfactory agreement, with negligible differences for all practical purposes. The model was equipped with four spherical Oqus markers, to ensure redundancy in case one marker was not detected by the cameras during the tests. We did not observe any failure in marker detection during the tests, which was as expected since the model was stationary.

The wave probes were calibrated four times during the testing period, where we found negligible differences in the calibration factor. The paired wave probes gave similar measurements when it was expected (for instance when the piston mode was excited).

The basin is equipped with a parabolic beach at the opposite side of the wavemaker. For the range of periods and wave steepnesses tested in the present study, some reflection is expected from the beach (typically 5% of the incident wave is expected to be reflected by the beach).

Reflection from the side walls were small, in general. Re-reflections from the wavemaker were observed, especially in beam sea conditions. This was accounted for during post-processing by using a time-window prior to re-reflection from the wavemaker and the parabolic beach (in regular waves at least).

Repetition tests were carried out in both irregular and regular waves for the configurations MP2 and MP3. The same irregular wave was repeated five times for two different wave steepnesses. The standard deviation is less than 1% relative to the mean value in the irregular wave tests. At narrow peaks close to the first sloshing mode for MP2 in beam sea the standard deviation of the free-surface motion in the moonpool is approximately 4% relative to the mean value. In regular waves, the same trend is observed, where the standard deviation in the repetition tests is less than 1% relative to the mean value, in general, but approximately 4.5% at the first sloshing mode with MP2. We believe this is related to the fact that swirling occurs in beam sea for MP2, which is discussed in § 4.2. In regular waves, the same wave was repeated eight times for selected wave periods.

Videos of the free-surface motion in the moonpool was recorded from two different angles. In addition, side-and under-water view video recordings of the model were also made. All of the video recordings were saved for each test automatically, where the video recordings started and ended at the same time as the other measurements. The video recordings served as a tool to interpret the results, such as swirling type sloshing inside the moonpool. We present snapshots and video recordings from the model tests herein.

Linear potential flow simulations

Numerical simulations assuming linear potential flow theory were carried out using the panel code WAMIT, which solves the radiation and diffraction problems and body motions in frequency domain. A lowest-order method is applied, i.e. the source strength over each panel is assumed to be constant. The effect of flow separation at sharp corners is not accounted for. Thus, we expect the resonant body and moonpool responses to be over-predicted where there is dominant viscous damping. However, we expect that the natural periods, and the responses outside resonance, to be well predicted.

Simulations

Meshes for the three geometries are created using a code developed for the present study. The code also generates the mesh of the internal free-surface, for removal of irregular frequencies [START_REF] Lee | Wamit user's manual, ver. 6.4[END_REF]. In addition to the geometries MP1, MP2 and MP3, an additional geometry without moonpool is included, referred to as MP0 (cf. Fig. 4). This model is included to discuss the effects caused by the moonpool on the vessel motions. Symmetry conditions are applied about the xz-plane at = y 0. The water depth is defined to be 5 m, i.e. the same as in the experiments.

Results for the free-surface elevation in the moonpool are extracted at the same locations as in the experiments using numerical wave probes.

The present simulations are carried out with approximately 5000 panels. The required CPU time for each wave period was approximately two minutes on a single 3.6 GHz CPU.

In the present simulations, we have applied empirical damping in roll to avoid unphysically large roll motions by following [START_REF] Ikeda | Roll damping of a sharp-cornered barge and roll control by a new-type stabilizer[END_REF]. Without the additional damping, clearly unrealistically large roll motions were predicted, also affecting sway motions.

Convergence study

Convergence studies were carried out for all four geometries. In general, the results converge faster for high wave periods, as expected. For the coarsest meshes and small periods, some unphysical peaks are observed in the RAOs. Satisfactory convergence is obtained for N ≳ 2000, where the difference relative to the finest mesh is less than 1%. In the present study, a mesh with approximately 5000 panels over the half body is adopted.

Results and discussions

Experimental and numerical amplitude-dependent RAOs of the ship and moonpool responses are presented in this section. We present the steady-state first harmonic response amplitudes made non-dimensional with respect to the first harmonic incident wave amplitude, except in § 4.4, where we present the higher harmonic responses. Moonpool responses as measured by WP1 (cf. Fig. 2) are presented in a body-fixed coordinate system. Experimental results in irregular waves with three different wave steepnesses and in regular waves with two different wave steepnesses are presented along with numerical results using WAMIT. The moonpool and heave responses are presented in head and beam seas. Pitch RAOs are presented in head sea. Surge and sway RAOs in head and beam seas, respectively, are presented for MP2 and MP3. Natural periods determined according to [START_REF] Molin | On the piston and sloshing modes in moonpools[END_REF] and from WAMIT simulations by solving the radiation problem are presented in Table 3, and will be used to discuss the results presented in this section. The piston mode periods are also calculated using theory presented by [START_REF] Molin | On natural modes in moonpools and gaps in finite depth[END_REF], which improves the estimates of the piston mode period, since the length and width of the barge are finite. Fig. 4. The panel models for the presently studied geometries using WAMIT. The figures present a coarse mesh (N ≈ 800) for the sake of clarity. The used mesh has N ≈ 5000. Symmetry conditions are applied about the xz-plane at = y 0. The mesh of the internal free-surface is not shown.

MP1

Fig. 5 presents moonpool, heave and pitch RAOs from experiments and numerical simulations for the MP1 configuration. In general, the agreement between experiments and numerical simulations is good outside resonance. Close to resonance, the response is over-predicted by WAMIT, as expected. The amplitude dependency is clear in the RAOs, which indicates that flow separation at sharp edges provides dominant damping on the piston mode response. The piston mode period predicted from the radiation problem in heave is presented in Table 3. Due to the hydrodynamic coupling between the vessel and the piston mode responses, the piston mode period will differ, but is = T 1.06 s in the experiment. This is very close to that found by solving the radiation problem. [START_REF] Fredriksen | Wave-induced response of a floating two-dimensional body with a moonpool[END_REF] found that this interaction effect was more significant, but in their study the

= V V / 0.25 MP S
. Here,

= V V / 0.015 MP S
for MP1. This indicates that the coupling between ship and moonpool responses is dependent on the size of the moonpool relative to the submerged ship volume. The moonpool response at = T 1.06 s is illustrated in Fig. 6.

Table 3

Natural periods according to [START_REF] Molin | On the piston and sloshing modes in moonpools[END_REF], and with WAMIT by solving the radiation problem and identifying the peaks in the RAOs. Natural periods for T < 0.55 s are not calculated (NC) with WAMIT. Around the piston mode, WAMIT significantly over-predicts the heave response in a narrow wave period range due to large moonpool response. The experiments suggest that the influence of the moonpool on the heave response is minor, with a small coupled resonance peak around the piston mode. The hydrodynamic interaction effect is not as pronounced as in [START_REF] Fredriksen | Wave-induced response of a floating two-dimensional body with a moonpool[END_REF] since the V MP /V S -ratio is much smaller, as discussed.

The moonpool's effect on the pitch motions is negligible for all practical purposes according to the experimental results. WAMIT predicts that this effect is more prominent in a narrow period range, since the piston mode response is over-predicted.

The experimental irregular and regular wave tests with = H/ 1/60 are compared in Fig. 5, where the RAOs from regular and irregular wave tests illustrate consistent results. This applies both for the vessel motions and the moonpool responses.

MP2

The moonpool and vessel motion RAOs for the MP2 configuration are presented in Fig. 7. In the proximity of the first sloshing mode, WAMIT greatly over-predicts the moonpool response, more so than for the piston mode. This indicates that the wave radiation damping is even smaller in the proximity of the first sloshing mode, than for the piston mode, as can be expected with the present ratio between the moonpool length and width, and ship draft.

The amplitude-dependent moonpool RAOs from the irregular wave tests in head sea illustrate an increase in the RAO with increasing wave steepness. This is contrary to what we expect, and the reason is that nonlinear effects are dominant. The RAOs from the regular wave tests in head sea in the vicinity of the first sloshing mode show significantly lower responses. This is discussed in more detail in Section 4.4.

In beam sea conditions, the RAOs from the irregular wave tests show similar trends to those observed in head sea. However, the RAOs in the proximity of the first sloshing mode from the regular wave tests are much more significant than in head sea. From the video recordings we observed that swirling occured in beam sea conditions as also confirmed by studying the phases between the WP1, WP2, WP5 and WP6 inside the moonpool (cf. Fig. 2). This is illustrated in Fig. 9. Swirling is a nonlinear effect, and is a combination of two antisymmetric standing waves [2, p. 4 and 135], which results in a progressive wave in the angular direction. In this study, we observe that swirling occurs in the counterclockwise direction. The direction of the angular motion is dependent on the transient conditions in the moonpool. We performed eight repetition tests at = T 0.7 s, and in each test we ob- served that swirling occured in the counterclockwise direction each time. Swirling induces loads on the ship and equipment in the moonpool both in the x-and y-directions, on contrast to a single standing wave, which only induces a force in the direction of the standing wave. Swirling in closed sloshing tanks might induce significant loads in the longitudinal and transverse direction. The present experimental results indicate that the influence of the swirling motion is small on the surge and sway motions, as illustrated in Fig. 8. This indicates that the loads induced in the longitudinal and transverse directions are too small to induce any significiant ship motions. However, the swirling motion might induce large loads on an object in the moonpool.

The heave response in head sea is clearly influenced by the moonpool response in the proximity of the piston mode, and vice versa. The hydrodynamic coupling between the heave and the moonpool motions results in significant influence of the moonpool on the heave RAOs in the proximity of the piston mode, much more so than for MP1. The reason is that MP2 configuration has four times larger V MP /V S -ratio. A small peak in the heave RAOs is also observed around the piston mode period in beam sea conditions (cf. Fig. 7d), although the model tests indicate that this peak is not as prominent as indicated by the WAMIT simulations.

Coupling between the moonpool and pitch responses is expected, since the center of the moonpool is located 0.16L from the center of the vessel. The pitch natural period is shifted towards a lower period by about 4% due to the coupled heave-pitch-piston mode motion, and the resonant pitch response is larger than without moonpool.

In Section 4.1, we discussed that the hydrodynamic interaction effect due to coupling between the ship and moonpool responses were small for MP1. However, for MP2 this coupling is significantly more pronounced. Especially in head sea conditions, where we observe significant influence of the moonpool in the heave RAOs. We also observe moderate influence on the pitch RAOs due to the piston mode oscillations. In full scale, the piston mode natural period corresponds to = T 6.8 s, which is in the re- levant range with respect to a real sea state. The predicted heave response from the experiments for MP2 is an order of magnitude larger than for a similar vessel without a moonpool in the proximity of the piston mode. This coupling effect can therefore influence the operability of the vessel, not only with respect to relatively large moonpool oscillations, but also significantly larger heave responses.

The regular and irregular wave tests illustrate consistent results for = H/ 1/60, as for MP1. In the proximity of the first sloshing mode, large differences are observed, which are discussed in more detail in Section 4.4. Fig. 9 illustrates the moonpool response in the proximity of the first transverse sloshing mode in beam sea for MP2.

MP3

The moonpool responses in head and beam sea are presented in Fig. 10a andb for the MP3 configuration. Since the moonpool is relatively large, several resonant modes are excited in the tested range of wave periods. The piston mode is excited at approximately 1.17 s, which is 10% lower than that predicted by WAMIT and [START_REF] Molin | On natural modes in moonpools and gaps in finite depth[END_REF] in forced heave conditions, and 22% lower than the prediction according to [START_REF] Molin | On the piston and sloshing modes in moonpools[END_REF] (cf. Table 3). The shift is, as mentioned earlier, due to the hydrodynamic coupling between the vessel and moonpool motions. The shift in the first longitudinal sloshing mode due to the coupling effect is negligible.

An interesting observation is that the piston mode period in freely floating conditions almost coincides with the first longitudinal sloshing mode period in freely floating conditions due to the significant coupling between ship and moonpool motions. The natural periods are so close that the two peaks in the moonpool RAOs are not separated (see ⓔ in Fig. 10a). We observe evidence of dominant damping due to flow separation at sharp edges at the moonpool inlet in the proximity of the piston and first longitudinal sloshing modes, where WAMIT greatly overpredicts the responses. The fact that the piston and first longitudinal sloshing modes almost coincide makes it difficult to distinguish which resonant moonpool mode contributes to the different hydrodynamic effects in the proximity of the piston and first longitudinal sloshing modes.

The second longitudinal sloshing mode is excited at approximately 1.04 s, with a small shift of 1% relative to the one found by solving the radiation problem (see ⓓ in Fig. 10a). The amplitude-dependency in the RAOs in the vicinity of the second longitudinal sloshing mode indicates that damping due to flow separation at sharp corners is dominant even for the second longitudinal sloshing mode. The wave radiation damping associated with the second longitudinal sloshing mode is small. This is observed both in head and beam sea conditions.

The third longitudinal sloshing mode is excited at approximately 0.89 s (see ⓒ in Fig. 10a), where it is clear that the irregular wave RAOs are larger than the regular wave RAOs in head sea conditions. The irregular wave RAOs increase with increasing wave steepness, which indicates that nonlinear effects are dominant in the proximity of the third longitudinal sloshing mode. The reason is secondary resonance, which excites the third longitudinal sloshing mode (cf. Section 4.4). The significant over-predition by WAMIT indicates that wave radiation damping is small, and that damping due to flow separation at sharp edges is dominant, in addition to the clear nonlinear free-surface effects.

The influence of the moonpool responses in the heave RAOs is significant. Fig. 10c andd indicate that both the piston and the second longitudinal sloshing modes affect the heave RAOs in both head and beam sea conditions considerably. Although the experimental peaks are not as pronounced as indicated by WAMIT, we observe substantial influence of the moonpool on the heave RAOs when compared to those for a ship without a moonpool. In particular, there is a pronounced heave resonance peak at the piston/first sloshing mode period, something which is hardly present for a vessel without moonpool. In the proximity of the second longitudinal sloshing mode, the heave RAOs are approximately six times larger than the one for a ship without a large moonpool. This will influence the operability of a vessel, more so than for MP1 and MP2. Similarly, in the proximity of the piston mode, the heave responses can be as much as three times larger than the heave response for MP0. The pitch RAOs in head sea conditions are presented in Fig. 10e. There is overall moderate influence of the moonpool on the pitch RAOs, except in the proximity of the piston and first longitudinal sloshing modes, where there is a coupled pitch resonance peak with 3-4 times the pitch motion compared to MP0. The moonpool's effect on the pitch RAOs is not as prominent as for the heave RAOs, which we believe is because the center of the moonpool is located close to the center of the ship.

The regular and irregular wave tests illustrate consistent results for = H/ 1/60, both for vessel and moonpool responses. When the moonpool response is dominated by nonlinear effects, some differences are observed, which is discussed in detail in §4. [START_REF] Guo | Resonant water motions within a recessing type moonpool in a drilling vessel[END_REF].

The first transverse sloshing mode is excited at = T 0.7 s, as illu- strated in Fig. 10b. The moonpool response at = T 0.7 s is significant, where the moonpool response is approximately six times the incident wave amplitude for the lowest wave steepness in regular waves. From Fig. 11 b we notice a peak in the proximity of the first transverse sloshing mode in the sway RAOs, which is due to the resonant fluid motion in the moonpool. The sway response is approximately five times larger with the large moonpool, relative to the response for a ship without a moonpool. This can be important to consider for a ship in operational conditions, although the peak in the sway RAOs is not as high and wide as in surge.

Fig. 11 a presents the surge RAOs in head sea conditions, where the RAOs suggest that the influence of the moonpool in the proximity of the piston and first longitudinal sloshing modes is significant. In this range, the surge RAOs are approximately five times larger than the RAO for the ship without moonpool (i.e. MP0). In addition to the first longitudinal sloshing mode, the piston mode influences the surge motion since the piston mode shape is nonuniform. If the ship is performing an operation that is sensitive to the surge motion of the vessel, the coupling effect between the moonpool and surge motions could have significant influence on the operability. The WAMIT simulations indicate large coupling between the narrow-banded resonant peaks and the surge responses in narrow wave period bands at = T 0.7 s and = T 0.89 s in head seas. However, the model tests suggest that this effect is negligible.

Assuming that such a large moonpool is built in reality, it could be imagined it would be with a purpose of lowering relatively large marine structures into the sea through the moonpool. In addition to the significant vertical loads that the object will be exposed to, large horizontal loads on the structure due to sloshing motion in the moonpool must also be considered. For instance, in a scenario where the first sloshing mode is excited, the largest vertical water motion and acceleration of the moonpool will be at the moonpool ends, while the largest horizontal motion and acceleration will be close to the center of the moonpool. This means that the large moonpool motion can induce large vertical and/or horizontal loads on the structure depending on where the structure is lowered through. For small-volume bodies the loads may to a first approximation be estimated by the Morison model with the (undisturbed) coupled moonpool response as ambient flow.

When discussing mode shapes, the term is generally associated with those found by solving the eigenvalue problem, such as in [START_REF] Molin | On the piston and sloshing modes in moonpools[END_REF]. The natural modes and corresponding mode shapes are then independent of the excitation. However, we discuss the difference in the piston mode shape due to the variation of the incident wave along the moonpool length (cf. Fig. 12). This is of relevance with respect to ship motions and loads on an object in the moonpool. Fig. 13 presents the piston mode shapes in forced heave, fixed ship with incident wave and freely floating conditions predicted by WAMIT for the MP3 configuration. The first harmonic amplitudes from the model tests in regular waves with = H/ 1/30 are presented along with the WAMIT simulations. The difference in the mode shape in fixed and forced heave conditions is negligible for MP1 and MP2. However, a significant difference is observed in the mode shape for MP3. In forced heave, the mode shape predicted by WAMIT is in good agreement with theory presented in [START_REF] Molin | On the piston and sloshing modes in moonpools[END_REF]. In fixed conditions, the piston mode shape changes significantly due to the variation in the excitation over the moonpool length due to the incident wave. This effect is negligible for MP1 and MP2, since the moonpool length, L m , is small compared to the wavelength, λ. The L m /λ-ratio in the proximity of the piston mode is approximately unity for MP3 ( = L / 0.97 m ), which gives a significant change in the piston mode shape for MP3 in waves. Traditionally, we use the term piston mode to refer to the zeroth mode in a moonpool. The term indicates that the moonpool oscillates nearly uniformly. The present results indicate that this is a good approximation for small moonpools, but is not so for a ship with a large moonpool in freely floating conditions.

The differences in the hydrodynamic interaction between ship and moonpool responses for MP1, MP2 and MP3 is exemplified in Fig. 14, where the heave RAOs from the irregular wave tests with wave steepness = H/ 1/80 are presented.

Higher harmonic moonpool RAOs

Selected higher harmonic amplitude-dependent RAOs of the moonpool response from the regular wave tests are presented in Fig. 15. The higher harmonic vessel motions are negligible for all purposes, and not presented herein. The nth harmonic moonpool . The third harmonic moonpool responses are larger than (kζ a ) 3 . However, the contribution of the third harmonic moonpool responses on the total response is small, whereas the second harmonic moonpool responses are large. We present amplitude-dependent moonpool RAOs of the second and third harmonic components. The harmonic components beyond = n 3 were found to be small and are not presented.

The second harmonic moonpool responses for MP2 in head sea are presented in Fig. 15a. Large second harmonic responses are observed in the proximity of ⓐ and ⓑ. The activation of a higher mode is called secondary resonance [2], and is a nonlinear effect where the primary mode excites a higher mode. The significant second harmonic moonpool responses for MP2 in head sea explain the nonlinear effects observed in the RAOs from the irregular wave tests in the proximity of = T 0.7 s (cf. Fig. 7a), where the RAOs clearly increase with increasing wave steepness. The first harmonic moonpool response in head sea in the proximity of the first longitudinal sloshing mode is relatively small from the regular wave tests (cf. Fig. 7a). The second harmonic moonpool response is, however, significant at = T 1.4 s. This is illustrated at ⓑ in Fig. 15a. Due to activation of secondary resonance in the proximity of = T 1.4 s, the first longitudinal sloshing mode is excited, which results in larger responses than the first harmonic moonpool responses at = T 0.7 s (cf. Fig. 16). Similar results were obtained for MP2 in beam sea.

The second harmonic moonpool responses for MP3 in head sea are presented in Fig. 15b. The second harmonic responses are significant. For instance, due to secondary resonance, the second longitudinal sloshing mode is excited in a regular wave with wave period = T 2.06 s. The third longitudinal sloshing mode is even more pronounced (cf. ⓓ in Fig. 15b).

The second harmonic moonpool responses for MP3 in beam sea are presented in Fig. 15c (cf. Fig. 17,). The first transverse sloshing mode is excited due to secondary resonance at = T 1.4 s (ⓒ in Fig. 15c). The third harmonic moonpool responses for MP3 in head sea are presented in Fig. 15d. The third harmonic moonpool responses are significant. The fifth longitudinal sloshing mode is excited due to wave excitation at = T 2.06 s. The third harmonic moonpool responses for MP2 are not as pronounced as for MP3 in the tested range of wave periods.

The amplitude-dependence in the higher harmonic moonpool RAOs in Fig. 15 can be attributed to several nonlinear effects. We believe the main candidate to the amplitude-dependence is dominant damping due to flow separation. For instance, at ⓓ and ⓔ in Fig. 15b, clear amplitudedependency is illustrated.

In a realistic scenario, a ship will operate in irregular waves. In such conditions, we observe that it is important to account for the fact that higher modes might be excited when the moonpool dimensions are large enough. Even for MP2, which has a moderate moonpool size, energy is transferred to higher modes. MP2 has moonpool dimensions 13.8 m × 13.8 m in full scale, which means that this is relevant with respect to existing ships. In full scale, the first longitudinal and transverse sloshing modes are approximately 4 s for MP2. In the North Sea, the waves will not contain significant energy in this range. However, the first longitudinal and transverse sloshing mode can be excited if there is enough energy in the wave spectrum around 8 s due to excitation of higher modes, which is in the range of periods with significant wave energy in a realistic sea environment. Linear theory will not capture such effects, since they are nonlinear.

Conclusions

We investigated the hydrodynamic interaction between moonpool and vessel responses for a realistic freely floating three-dimensional vessel. Dedicated model tests were carried out in the Ocean Basin at Sintef Ocean in Trondheim, Norway. The ship and moonpool responses were investigated in head and beam sea conditions in a relevant range of periods for a vessel in the North Sea. The experiments were carried out in both regular and irregular waves. Numerical simulations assuming linear potential flow theory were carried out using WAMIT. Four different configurations were simulated, referred to as MP0, MP1, MP2 and MP3 (MP0 was without a moonpool and was not tested in the experiments). The V MP /V S ratios were 0.015, 0.064 and 0.43 for MP1, MP2 and MP3, respectively, where V S and V MP are the submerged ship and moonpool volumes. Resonant moonpool responses, and their influence on the ship motions, were discussed in terms of amplitude-dependent RAOs.

The hydrodynamic interaction between ship and moonpool responses for MP1 was negligible, while the influence of the moonpool on the vessel motions was significant for MP2, in particular, for the heave response. For MP3, the surge, sway, heave and pitch responses were strongly influenced by the moonpool. For instance, the heave motion was approximately ten times larger than for a ship without a moonpool in the proximity of the piston mode. Similarly, an increase in the heave RAOs of approximately six times was observed in the proximity of the second sloshing mode due to the coupling between ship and moonpool responses. Several sloshing modes were excited in the tested range of wave periods for MP3. The first transverse sloshing mode affected the sway motion in beam sea, where the experimental RAOs were approximately five times larger than for MP0. The piston and the first longitudinal sloshing modes affected the surge motion considerably. The piston and first longitudinal sloshing modes were not distinguishable.

In beam sea conditions with MP2, swirling occured in the proximity of the first transverse sloshing mode. Swirling is a nonlinear effect. Although the swirling motion was significant, it did not induce prominent surge and sway motions. However, swirling can be of concern with respect to operability of the moonpool.

The model tests suggested that secondary resonance was important to consider for the MP2 and MP3 configurations, since the higher modes were excited in a wave period range where a realistic sea environment contains significant wave energy. The first harmonic response in regular waves with MP2 at the first longitudinal sloshing mode in head seas was relatively small. The second harmonic response was, however, significant. We believe the main candidate to the amplitude-dependence in the first and second harmonic moonpool RAOs is the dominant damping due to flow separation at the moonpool entrance for both harmonics. In contrast to the traditional moonpool case where damping due to flow separation leads to lower amplitude-dependent RAOs for increasing wave steepness, amplitude-dependent RAOs from the irregular wave tests increase with increasing wave steepness due to secondary resonance. This effect is increasingly important when the moonpool dimensions increase.

Fig. 1 .

 1 Fig. 1. The ship model in the present experiments in the Ocean Basin. The model is 4 m long, the mid-ship beam is 0.8 m midships and the total height is 0.6 m. Different moonpool sizes can be tested by substituting boxes in the large opening.

Fig. 2 .

 2 Fig. 2. Bird's eye view of the wave probes (WP) in the moonpool for the present configurations. The model dimensions are also indicated.

Fig. 3 .

 3 Fig. 3. Measured and desired Pink Noise for Water Waves spectra with three different wave steepnesses.

Fig. 5 .

 5 Fig. 5. Responses in head and beam sea for MP1. Moonpool responses at the location of WP1 in a body-fixed coordinate system. Results for MP0 (without moonpool) are presented as pink dashed lines. The numbers shown near the peaks in (a-c) indicate the peak values from the WAMIT simulations.

Fig. 5 .

 5 Fig. 5. (continued)

Fig. 5 .Fig. 6 .

 56 Fig. 5. (continued)

Fig. 7 .

 7 Fig. 7. Selected responses in head and beam sea for MP2. Moonpool responses at the location of WP1 in a body-fixed coordinate system. Results for MP0 (without moonpool) are presented as pink dashed lines.

Fig. 7 .Fig. 7 .

 77 Fig. 7. (continued)

Fig. 8 .

 8 Fig. 8. Surge and sway RAOs in head and beam sea for MP2. Results for MP0 (without moonpool) are presented as pink dashed lines.

Fig. 9 .

 9 Fig. 9. Four snapshots from video recordings of the moonpool response for MP2 in beam sea, where swirling occured at = T 0.7 s in regular waves. See Video 2 for full video of the moonpool oscillations.

Fig. 10 .

 10 Fig. 10. Responses in head and beam sea for MP3. Moonpool responses at the location of WP1 in a body-fixed coordinate system. Results for MP0 (without moonpool) are presented as pink dashed lines.

Fig. 10 . 1 .

 101 Fig. 10. (continued)

Fig. 10 .Fig. 11 .

 1011 Fig. 10. (continued)

Fig. 12 .

 12 Fig. 12. Two snapshots from video recordings of the piston mode response for MP3 in head sea in regular waves at = T 1.15 s. See Video 3 for full video of the moonpool oscillations. The time instants are approximately T/2 apart, where T is the wave period. The arrows point to the maximum piston mode response in the moonpool. .

Fig. 13 .Fig. 14 .

 1314 Fig. 13. Piston mode shapes as predicted by WAMIT in forced heave, fixed ship in incident waves and freely floating conditions. The mode shapes are for MP3 in head seas at = y 0, where the mode shapes from WAMIT are established using 1000 evenly distributed points in the moonpool. The experimental mode shapes are from the regular wave tests with H/λ = 1/30, where we present the first harmonic response amplitudes. The mode shapes are made non-dimensional with respect to the maximum response in the moonpool, ζ max .Experimental and numerical results in freely floating conditions at = T 1.17 s are presented. In ship and forced heave, the wave and forced heave periods are both = T 1.30 s.

Fig. 15 .

 15 Fig. 15. Higher harmonic moonpool RAOs at WP1 in head and beam sea conditions. The RAOs are presented in a body-fixed coordinate system. The circles indicate which sloshing mode is excited due to secondary resonance.

Fig. 15 .Fig. 16 .

 1516 Fig. 15. (continued)

Fig. 17 .

 17 Fig. 17. Two snapshots from video recordings of the moonpool response for MP3 in head sea in regular waves at = T 1.77 s. Dominant excitation of the third longitudinal sloshing mode is observed. See Video 5 for full video of the moonpool oscillations. The time instants are approximately T/2 apart, where T is the wave period.

Table 1

 1 Main properties of the present ship model. A model scale of 1:34.5 is imagined. The moonpool center is measured relative to the origin of the coordinate system.

		Model scale	Full scale
	Length (L)	4.0 m	138 m
	Beam (B)	0.8 m	27.6 m
	Draft (D)	0.2 m	6.9 m
	Dimensions with MP1 (V MP /V S = 1.5%)
	Moonpool length (L m )	0.2 m	6.9 m
	Moonpool width (B m )	0.2 m	6.9 m
	Moonpool center ((c x , c y , c z ))	(0.65 m, 0 m, 0 m)	(22.4 m, 0 m, 0 m)
	Ship mass (m)	526.7 kg	21627.8 t
	Center of gravity (COG)	(-0.11 m, 0 m,	(-3.8 m, 0 m, 1.9 m)
		0.055 m)
	Radius of gyration ((r x , r y , r z ))	(0.28 m, 1.02 m,	(9.7 m, 35.2 m, 26.2 m)
		0.76 m)
	Dimensions with MP2 (V MP /V S = 6.4%)
	Moonpool length (L m )	0.4 m	13.8 m
	Moonpool width (B m )	0.4 m	13.8 m
	Moonpool center ((c x , c y , c z ))	(0.65 m, 0 m, 0 m)	(22.4 m, 0 m, 0 m)
	Ship mass (m)	502.69 kg	20642.3 t
	Center of gravity (COG)	(-0.15 m, 0 m,	(-5.2 m, 0 m, 2.0 m)
		0.057 m)
	Radius of gyration ((r x , r y , r z ))	(0.29 m, 1.024 m,	(10.0 m, 35.2 m,
		0.75 m)	25.9 m)
	Dimensions with MP3 (V MP /V S = 42.7%)
	Moonpool length (L m )	2.0 m	69 m
	Moonpool width (B m )	0.4 m	13.8 m
	Moonpool center ((c x , c y , c z ))	(-0.15 m, 0 m, 0 m)	(-5.2 m, 0 m, 0 m)
	Ship mass (m)	374.69 kg	15386.1 t
	Center of gravity (COG)	(-0.08 m, 0 m,	(-2.76 m, 0 m, 2.9 m)
		0.083 m)
	Radius of gyration ((r x , r y , r z ))	(0.30 m, 1.11 m,	(10.4 m, 38.3 m,
		1.13 m)	39.0 m)

Table 2

 2 Tested irregular pink noise waves in the ocean basin. Each irregular wave test was run for 15 minutes. The significant wave height, H s , is also presented.

	Spectrum type	H/λ	H s specified	H s achieved
	Pink noise for water waves	1/100	0.034 m	0.035 m
	Pink noise for water waves	1/80	0.046 m	0.049 m
	Pink noise for water waves	1/60	0.092 m	0.103 m

Declaration of Competing Interest

None.

Supplementary material

Supplementary material associated with this article can be found, in the online version, at 10.1016/j.apor.2019.102010


Daniel De Wolf
email: daniel.dewolf@univ-littoral.fr

Yves Smeers
email: yves.smeers@uclouvain.be

Generalized derivatives of the optimal value of a linear program with respect to matrix coefficients

Keywords:

 to the cases where derivatives may not be defined because of the existence of multiple primal or dual solutions.

Introduction

In the framework of linear programming, we consider the problem of estimating the variation of the objective function resulting from changes in some matrix coefficients. Our objective is to extend results already available for the right-hand side to this more general problem.

The interpretation of the dual variables as derivatives of the optimal value of the objective function with respect to the elements of the right-hand side is well known in mathematical programming. This result can be extended to the case of multiple dual solutions. The set of all dual solutions is then the subdifferential of the optimal value of the objective function, seen as a convex function of the right-hand side. The object of this paper is to extend these well known results to the derivative of the optimal value of the objective function with respect to matrix coefficients.

It is easy to show on a simple example that the objective function value of a linear program is not a convex function of the matrix coefficients. The subdifferential concept is thus inappropriate here. One must therefore resort to Clarke's notion of a generalized gradient. A characterization of this generalized gradient will be derived and sufficient conditions of existence of the generalized gradient will be given for this particular application of nonsmooth analysis.

The paper is organized as follows. Section 2 presents a study of the literature about generalized derivatives and subdifferentials of optimal function value. Then Section 3 presents the basic definitions and mains properties of nonmooth analysis that will be useful for our application. Then In Section 4 we recall the result of [START_REF] Freund | Postoptimal analysis of a linear program under simultaneous changes in matrix coefficients[END_REF] in the smooth case, namely the gradient of the optimal value function in the case where the optimal primal and dual solutions of the linear problem are unique. Then in Section 5 , a complete characterization of the generalized gradient for the case where the primal or dual solutions are not unique is established. Since the Locally Lipschitz property plays an essential role in this characterization, we will give sufficient conditions to prove the Local Lipschitz property. Section 6 presents a practical application of this characterization coming from the gas industry, namely, the two-stage problem of optimal dimensioning and operating of a gas transmission network. Then Section 7 presents some conclusions.

Generalized derivatives and subdifferentials of optimal value functions

Before establishing the main result concerning the generalized gradient of the optimal value function of a linear program with respect to matrix coefficients, let us say a word about the literature on the problem of the computation of the generalized derivatives and subdifferentials of the optimal value function in general optimization problems.

One of the first papers on the subject was the paper of [START_REF] Gauvin | The generalized gradient of a marginal function in mathematical programming[END_REF] who considered a general mathematical programming problem with equality and inequality constraints and pertubation of the right-hand side of the contraints, noted u i for constraint i . He estimated the generalized gradient of the optimal value of the

z(u ) = max x f (x ) , x ∈ R n s.t. g i (x ) ≤ u i , i = 1 , . . . m h i (x ) = u i , i = m + 1 , . . . p
Note that this case corresponds, for the particular case of a linear program, to a perturbation of the right-hand side b of the linear programming problem since we consider a perturbation (noted A in our case) that appears as a coefficient of the variables x of the problem .

Several developments concerning this case (i.e. the perturbation only in the right-hand side of the constraints) where made since that time. For example, the regularity properties of the optimal value function in nonlinear programs where the perturbation parameters appears only in the right-hand side was done by [START_REF] Craven | Regularity properties of the optimal value function in non linear programming[END_REF] . For affine case, namely when g(x ) = Ax + b, an expression is given for its directional derivative, not assuming the optimum to be unique. Recently, Höffner, Khan, and [START_REF] Höffner | Generalized derivatives of dynamic systems with a linear program embedded[END_REF] consider the computation of generalized derivatives of dynamic systems with a linear program embedded. They consider the optimal value of a linear program as a function of the right-hand side of the constraints and present an approach to compute an element of the generalized gradient. The approach is illustrated through a large-scale dynamic flux balance analysis example. Last year, [START_REF] Gomez | Generalized derivatives of lexicographic linear programs[END_REF] studied the generalized derivatives problem of parametric Lexicographic Linear programs using the lexicographic directional derivative (See [START_REF] Barton | Computationally relevant generalized derivatives: theory, evaluation and applications[END_REF] for a survey on the lexicographic directional derivative). The paper derives generalized derivatives information.

Unfortunately, this results can not be applied to our problem since the perturbation appears only in the right-hand side of the contraints and doesn't appear as a coefficient of the decision's variables. [START_REF] Rockafellar | Directional differentiability of the optimal value function in a nonlinear programming problem[END_REF] considered the directional derivative of the optimal value function in nonlinear programming problem with a perturbation u that appears in the left-hand side of the constraints:

z(u ) = max x f (x, u ) , x ∈ R n s.t. g i (x, u ) ≤ 0 , i = 1 , . . . m
(1)

h i (x, u ) = 0 , i = m + 1 , . . . p
Under the assumption that every optimal solution x * satisfies the second order constraints qualification condition , Rockafellar proved that the function z ( u ) is locally Lipchitz and finite. Rockafellar also gives an upper bound on the generalized derivative of Clarke of the function z ( u ). Note also that many developments were made from this original paper for the general case where the perturbation appears in the left-hand side of the constraints. For example, [START_REF] Thibault | On subdifferentials of optimal value functions[END_REF] considered a general mathematical programming problem in which the constraints are defined by multifunctions and depend on a parameter u . A special study is done of problems in which the multifunctions defining the constraints take convex values. For these problems, generalized gradients of z ( u ) are given in terms of the generalized gradients of the support functions of the multifunctions. Bonnans and Shapiro (20 0 0) studied the first order differentiability analysis of the optimal value function as a function of a parameter that appears in the objective function and in the lefthand side of the contraints. Under a constraint qualification condition, they give an upper bound on the directional derivative. Note that in our case of a linear program, we will give a complete characterization of the generalized gradient, not only upper bound on the directional derivative. [START_REF] Penot | Differentiability properties of optimal value functions[END_REF] considers the differentia-bility properties of optimal value functions for the particular case where the perturbation parameter only appears in the objective function. More recently, [START_REF] Mordukhovich | Subgradients of marginal functions in parametric mathematical programming[END_REF] consider the subgradient of marginal functions in parametric mathematical programming. The authors show that the subdifferential obtained for the corresponding marginal value function are given in terms of Lagrange multipliers. Last year, [START_REF] Im | Sensitivity analysis and robust optimization: A geometric approach for the special case of linear optimization[END_REF] studied the sensitivity analysis for the special case of linear optimization. In particular, he gives conditions for the objective function value of a linear problem to be a Locally Lipschitz function of matrix coefficients. We will use these conditions in our main characterization of the generalize gradient.

In the present paper, we shall give a complete characterization of the generalized gradient for a particular case, namely the linear case, and not only upper bound on directional derivative.

Basic definitions and properties

This section recalls some basic concepts and properties of nonsmooth optimization useful for our application. An introduction to the first-order generalized derivative can be found in [START_REF] Clarke | Optimization and nonsmooth analysis[END_REF] for the case of a locally Lipschitz function.

Definition 3.1. A function f from R n (or a subset of R n ) into R is locally Lipschitz if for any bounded set B from the interior of the domain of f there exists a positive scalar K such that | f (x ) -f (y ) | ≤ K x -y ∀ x, y ∈ B
where |.| denotes the absolute value and . the usual Euclidian norm.

The locally Lipschitz property can be interpreted as a finite bound on the variation of the function. It is well known that the locally Lipschitz property implies the continuity of f .

The Rademacher theorem says that a locally Lipschitz function f has a gradient almost everywhere (i.e. everywhere except on a set Z f of zero (Lebesque) measure on R n ). Definition 3.2. In the locally Lipschitz case, the generalized gradient is defined as the convex hull of all the points lim ∇f ( x k ) where { x k } is any sequence which converges to x while avoiding the points where ∇f ( x ) does not exist:

∂ f (x ) = con v lim k →∞ ∇ f (x k ) : x k → x, ∇ f (x k ) exists (2)
where con v denotes the convex hull.

Another essential concept in nonsmooth optimization is the directional derivative. This notion can also be generalized to the nonconvex case.

Definition 3.3. The generalized directional derivative of f evaluated at x in the direction d is defined (using the notation of Clarke) as

f 0 (x ; d) = lim sup y → x ;t↓ 0 f (y + td) -f (y ) t
In the convex case, this notion reduces to the classical notion of directional derivative

f (x ; d) = lim t↓ 0 f (x + td) -f (x ) t
We shall also use the following proposition for the proof of our characterization of the generalized gradient: Proposition 3.1. Let f be a function from R n into R almost everywhere continuously differentiable. Then f is continuously differentiable at x if and only if ∂f ( x ) reduces to a singleton.

Proof. See [START_REF] Clarke | Optimization and nonsmooth analysis[END_REF] .

We have the following general remark. The definition of the generalized gradient (2) is only valid for the locally Lipschitz case. If the function is simply almost everywhere differentiable, one can construct examples for which the generalized gradient is not defined. A more general definition based on the cone of normals is given by [START_REF] Clarke | Optimization and nonsmooth analysis[END_REF] in the case of a lower semi-continuous function.

Gradient of the optimal value function

Returning now to our problem, we consider the optimal value of a linear problem as a function of the matrix coefficients:

z(A ) = max x c T x subject to Ax = b x ≥ 0 ( 3
)
where c is the n -column vector of objective coefficients, x is the n -column vector of variables, A is the m × n matrix of left-hand side coefficients and b is the m -column vector of right-hand side coefficients.

We first recall the result for the smooth case. It is established in [START_REF] Freund | Postoptimal analysis of a linear program under simultaneous changes in matrix coefficients[END_REF] under the two following assumptions:

(H1) The optimal solution of the primal problem (3) is unique; (H2) The optimal solution of the dual problem of (3) is unique.

The result in the smooth case can then be written as follows:

Proposition 4.1. [START_REF] Freund | Postoptimal analysis of a linear program under simultaneous changes in matrix coefficients[END_REF] ). If assumptions (H1) and (H2) are both satisfied for A , then z ( A ) is continuously differentiable and we have that

∂z ∂a i j = -u * i x * j ( 4
)
where u * i is the optimal dual variable associated to row i and x * j is the optimal primal variable associated to column j.

A complete analysis of the subject in the differentiable case can be found in [START_REF] Gal | Postoptimal analyses[END_REF] .

Generalized gradient characterization.

Before examining the case where the optimal basis is not unique , we show on an example that z ( A ) does not enjoy any convexity property. Consider the following linear problem with a single parametric matrix coefficient:

z(a ) = min x x 1 + x 2 s.t. ax 1 + x 2 = 1 x 1 , x 2 ≥ 0
Using the constraint, x 2 can be substituted:

z(a ) = min x 1 + (1 -a ) x 1 s.t. 0 ≤ x 1 ≤ 1 a
The optimal objective function can thus be written explicitly as:

z(a ) = 1 if a < 1 1 /a if a ≥ 1
It is clear that z ( a ) is neither convex nor concave. Because of this lack of convexity, the notion to be used is the Clarke's generalized gradient.

If A is such that the linear program is infeasible, we define

z(A ) = -∞ .
Denote by dom ( z ), the domain where z ( A ) is finite. Before stating the characterization of the generalized gradient, we first recall the following propositions which result from [START_REF] Renegar | Some perturbation theory for linear programming[END_REF] .

Proposition 5.1. If the set of optimal primal solutions for A is unbounded, then A is not an interior point of dom ( z ) .

Proposition 5.2. If the set of optimal dual solutions for A is unbounded, then A is not an interior point of dom ( z ) .

We will use the following notation u × x for the outer product of an n -column vector u by the n -row vector x T . The following theorem states a characterization of the generalized gradient.

Theorem 5.1. If A is an interior point of dom ( z ) and if z ( A ) is locally Lipschitz in a neighborhood of A , then ∂z(A ) = con v {-u × x where u is any optimal dual solution and x is any primal optimal solution of (3) }

Proof.

1. Suppose first that there is a single optimal basis. Since A is an interior point of dom ( z ), we know by Propositions 5.1 and 5.2 that there are no extreme rays of primal or dual optimal solutions. In this case, a single optimal basis is a sufficient condition to have primal and dual nondegeneracy. We know from Proposition 4.1 that ∂z ( A ) reduces to a single matrix, which can be computed by the following formula:

∂z(A ) = { -u × x }
where u and x are the dual and primal solutions associated to the unique optimal basis for (3) . This proves the theorem. 2. Suppose next that there are several optimal basis. Lets first introduce some useful notation from linear programming. For a particular optimal basis, we denote by B the columns of matrix A corresponding to the basic variables noted x B and we denote by N the columns of matrix A corresponding to non basic variables, noted x N . The constraints can by rewritten as follows:

Ax = b ⇔ (B, N)(x B , x N ) T = b ⇔ Bx B = b ⇔ x B = B -1 b
Denote c b the objective coefficient corresponding to basic variables, it follows that the objective function can be rewritten for a particular basis as follows:

z(A ) = c T B B -1 b
We first prove the following inclusion:

con v {-u × x such that u T = c T B B -1 , x B = B -1 b, x N = 0
and B corresponds to an optimal basis of (3) } ⊂ ∂z(A )

Let B corresponds to an optimal basis of (3) . Since this optimal basis is not unique, there must be at least one non-basic variable x j with a zero reduced cost:

c j -c T B B -1 a j = 0
where a j denotes the column j of matrix A . Using the definition of the vector of dual variables u T = c T B B -1 , this condition can be rewritten:

c j -u T a j = 0
We can exclude the pathological case where u = 0 . In fact, this case can be treated by a perturbation of the objective coefficients. This only requires to consider A as the extended matrix of the system where the objective is added as a constraint. We can thus take u i different from zero and define the following perturbation of the column a j for any column with zero reduced cost: if u i > 0, subtract > 0 ( < 0 if u i < 0) from the i th component of a j . For the perturbed problem, all the reduced costs are strictly positive, and therefore the optimal solution becomes unique.

More specifically, let a sequence of matrices A ( ) where the objective function is differentiable conver ging towards matrix A where it is not differentiable. The primal and dual optimality conditions hold for each point of the sequence and the optimal primal and dual variables are continuous functions of the perturbed matrix since this latter is a perturbation of an invertible matrix at the limit point. The limits of the primal and dual variables at each point of the sequence thus exist and satisfy the primal dual relations at the limit point. It can be concluded that all the points associated with all the optimal basis for the problem (3) belong to the generalized gradient, since they can be written as the limit of perturbed problems where the optimal basis is unique and so where the function z is continuously differentiable in A ( ) for all > 0. Now consider the inverse inclusion:

∂z(A ) ⊂ con v {-u × x such that u T = c T B B -1 , x B = B -1 b, x N = 0
and B corresponds to an optimal basis of (3) }

Let { A k } be any sequence such that ∇z ( A k ) exists and converges.

The optimal basis associated with each A k does not have to be unique. As shown by [START_REF] Freund | Postoptimal analysis of a linear program under simultaneous changes in matrix coefficients[END_REF] , we can have, for a given matrix A , several degenerate optimal basis although z is continuously differentiable for matrix A . We will show, in this case, that any optimal basis associated with A k must give the same point ( . . . , -u i x j , . . . ) .

Suppose the opposite, that is, there are two different optimal basis for A k , giving two different points ( . . . , -u 1 i x 1 j , . . . ) and ( . . . , -u 2 i x 2 j , . . . ) respectively. As done in part i) of the proof, the matrix A k can be perturbed in order to have no more than one of the two basis optimal. By taking this limit, we obtain that the first point is in the generalized gradient. By applying the same procedure of perturbation to the second basis, we show that the second point is also in the generalized gradient. We can therefore conclude that ∂z ( A k ) is not a singleton. Applying Proposition 3.1 , this contradicts the fact that z is continuously differentiable in A k . The gradient can therefore be associated with any of the optimal basis. Note by { β k } a sequence of optimal basis for A k (i.e. β k is an optimal basis for A k ). By a basis β, we mean here a partition of the variables between basic and non-basic variables. As { β k } is an infinite sequence of basis and as there is only a finite choice of m columns among the n columns of the A matrix, so there must be a special basis β which is repeated infinitely often in the sequence. Let { B l } be the subsequence corresponding to this basis which is repeated infinitely often.

The corresponding subsequence { ( . . . , -x j u i , . . . ) } l of gradients associated with this basis converges to the same point as the original sequence. As

c T N -c T B (B l ) -1 ≥ 0 (B l ) -1 b ≥ 0
for all l , these inequalities remain true for l → ∞ and so { B l } converges to an optimal basis for (3) . This completes the proof of the reverse inclusion. 3. We finally show that ∂z(A ) = con v {-u × x where u is any optimal dual solution and x is any primal optimal solution of (3) } Because the sets of primal and dual optimal solutions corresponding to point A are bounded by Propositions 5.1 and 5.2 , u and x are convex combinations of extreme dual and primal solutions respectively. Let

u = k μ k u k where k μ k = 1 and μ k ≥ 0 x = l λ l x l where l λ l = 1 and λ l ≥ 0
Suppose first that u is a convex combination of extreme u k while x is an extreme optimal point. One has

u i x j = k λ k u k i x j
for a given set of λ k and for all i and j . Therefore

-u × x = - k λ k u k × x k
This implies that con v {-u × x where u is any optimal dual solution and

x B = B -1 b, x N = 0 , where B is the optimal basis } = con v {-u × x where u = c T B B -1 , x B = B -1 b,
x N = 0 and B is any optimal basis of (3) }

The same reasoning can be made in order to relax the requirement that x is an extreme solution into the weaker one that x is any optimal solution of problem (3) .

Before illustrating the theorem on an example, let us say a few words about the requirements for A to be an interior point of dom ( z ) and for z ( A ) to be Lipschitz in a neighborhood of A . [START_REF] Im | Sensitivity analysis and robust optimization: A geometric approach for the special case of linear optimization[END_REF] proves that these two requirements holds true if the following conditions are satisfied: Assumption 5.1. The matrix A is of full rank and the Slater constraints qualification is satisfied.

Proposition 5.3. If matrix A is a full rank matrix and if the Slater constraints qualification is satisfied, then

• A is an interior point of dom ( z ) and • the function z ( A ) is Lipchitz in a neighborhood of A.

Proof. See [START_REF] Im | Sensitivity analysis and robust optimization: A geometric approach for the special case of linear optimization[END_REF] , pages 74-76.

As indicated by [START_REF] Höffner | Generalized derivatives of dynamic systems with a linear program embedded[END_REF] , any linear program can be reduced to an equivalent linear program that satisfies the full rank property for A by removing linearly dependent rows.

The following simple example illustrates Theorem 5.1 :

z(a ) = max x x 1 + x 2 s.t. x 1 + 2 x 2 ≤ 3 x 1 + ax 2 ≤ 2 x 1 , x 2 ≥ 0
The feasible region and the objective function are represented in Fig. 1 for the particular choice of a = 1. For a = 1 , there exists two different basic solutions. The first one is obtained with x 1 and x 2 in the basis: (x 1 , x 2 ) = (1 , 1) and the reduced cost of s 1 , the first slack variable, is zero. The second solution is obtained by taking x 1 and s 1 in the basis: (x 1 , x 2 ) = (2 , 0) and the reduced cost of x 2 is zero. In both cases, the optimal dual values are given by (u 1 , u 2 ) = (0 , 1) .

Take a = 1and let go to zero. We obtain the first solution and the reduced cost associated to s 1 is strictly negative. Take a = 1 + and let go to zero. We obtain the second solution and the reduced cost associated to x 2 is strictly negative. The extreme points of the generalized gradient are thus:

-u 2 x 2 = -1 (first case) -u 2 x 2 = 0 (second case).

One therefore obtains: In fact, the general expression of z ( a ) can be computed explicitly as:

∂z(1) = [ -1 , 0]
z(a ) = 3 a -5 a -2 if a < 1 2 if a ≥ 1
The graph of the optimal value of the function z ( a ) is represented in Fig. 2 as a function of parameter a . The two points -1 and 0 correspond thus to the left-and right-derivatives of z ( a ) at point a = 1 respectively.

Pratical application

The motivation for considerating the computation of the generalized gradient of the objective function of a linear problem with respect to the matrix coefficients is the general problem where there is a two-stage problem where at the first stage a capacity investment decision is made and at the second stage the operating of the system is optimized taking into account this investment decision. In some cases, the investment decision appears in the righthand side (such as capacity level decision). We consider the case where the investment decision appears as coefficient of the secondstage decision variables. Let us illustrate this fact with an example from the gas industry.

Consider the problem of the optimal dimensioning of pipe networks for the transmission of natural gas. See, for example, De Wolf and Smeers (1996) who consider the following two-stage problem for investment and exploitation of gas transmission networks. At the first stage, the optimal diameters of pipe lines , denoted D , must be determined in order to minimize the sum of the direct investment cost function , denoted C ( D ), and Q ( D ), the future operating cost function .

min D F (D ) = C(D ) + Q (D ) s.t. D i j ≥ 0 , ∀ (i, j) ∈ SA (5)
where SA denotes the set of arcs of the network.

The operations problem for a given choice of the diameters can thus be formulated as follows:

Q (D ) = min f,s,p j∈ N s c j s j s.t. ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ j | (i, j ) ∈ A f i j - j| ( j,i ) ∈ A f ji = s i ∀ i ∈ N sign ( f i j ) f 2 i j = K 2 i j D 5 i j (p 2 i -p 2 j ) ∀ (i, j) ∈ SA s i ≤ s i ≤ s i ∀ i ∈ N p i ≤ p i ≤ p i ∀ i ∈ N (6)
where the variables of the problem are f ij , the flow in the arc ( i , j ), s i , the net supply at the node i and p i , the pressures at node i . The set of nodes is denoted N . For simplicity of notation, we define the variable π i as the square of the pressure at node i :

π i = p 2 i .
Let us replace in the only nonlinear relation of exploitation problem (6) , the D ij variable by the following substitute:

x i j = D 5 i j

In fact, taking the x ij as parameters, we find that they appear as linear coefficients of the squared pressure variables in the equation:

sign ( f i j ) f 2 i j -K 2 i j D 5 i j (π i -π j ) = 0 .

(7)

De Wolf and Smeers (20 0 0) solve the problem of the gas transmission problem by an extension of the Simplex algorithm using piecewise linearisation of the first term sign ( f i j ) f 2 i j . We are thus back in the linear case. Let w * i j be the optimal value of the dual variable associated with constraint (7) . Applying Theorem 5.1 , one obtains an element of the generalized gradient by:

∂Q ∂x i j = w * i j K 2 i j (π * i -π * j ) (8)
Now, to obtain an element of the generalized gradient with respect to the original variables ( D ij ), one uses the chain rule for the composition of derivative with:

∂x i j ∂D i j = 5 D 4 i j .
It is then easy to prove that the following expression gives an element of the generalized gradient of Q ( D ):

∂Q(D ) ∂D i j = w * i j (π * i -π * j )5 K 2 i j D 4 i j (9)
This formula gives thus an element of the generalized gradient, which is the only information required by the bundle method (See [START_REF] Lemaréchal | Nondifferentiable optimization[END_REF] ) used to solve the two stage problem. See De [START_REF] De Wolf | Optimal dimensioning of pipe networks with application to gas transmission networks[END_REF] for the application of the bundle method to this two-stage problem.

Conclusions

It has been shown in this paper how the first-order derivatives of the optimal solution of a linear program with respect to matrix coefficients can be generalized to the nonsmooth case, even when the optimal function as a function of matrix coefficients admits breakpoints. Our result, Theorem 5.1 , emphasizes the fundamental role played by bases in this respect. The extreme points of the generalized gradient correspond to all the different optimal basis. A practical application to gas transmission network optimization, which was in fact the motivation for considering such a formula, was then presented.

  https://doi.org/10.1016/j.ejor.2019.11.020 0377-2217/© 2020 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license ( http://creativecommons.org/licenses/by-nc-nd/4.0/ ) problem considered as a function of the right-hand side perturbation:

Fig. 1 .

 1 Fig. 1. Illustrative example.

Fig. 2 .

 2 Fig. 2. Graph of the optimal value of the objective function.

Acknowledgments

We would like to thank Professor Guy de Ghellinck for valuable discussions and suggestions. We would also like to thank Professors Jean-Jacques Strodiot and Claude Lemaréchal for insightful discussions on the state of the art in nonsmooth optimization. Finally, we would like to thank the two anonymous reviewers for valuable suggestions to improve the paper.