# Adaptive Model Predictive Control for Engine-Driven Ducted Fan Lift Systems
using an Associated Linear Parameter Varying Model

Hanjie Jiang  Ye Zhou<EMAIL_ADDRESS>Hann Woei Ho  Wenjie Hu

###### Abstract

Ducted fan lift systems (DFLSs) powered by two-stroke aviation piston engines
present a challenging control problem due to their complex multivariable
dynamics. Current controllers for these systems typically rely on
proportional-integral algorithms combined with data tables, which rely on
accurate models and are not adaptive to handle time-varying dynamics or system
uncertainties. This paper proposes a novel adaptive model predictive control
(AMPC) strategy with an associated linear parameter varying (LPV) model for
controlling the engine-driven DFLS. This LPV model is derived from a global
network model, which is trained off-line with data obtained from a general
mean value engine model for two-stroke aviation engines. Different network
models, including multi-layer perceptron, Elman, and radial basis function
(RBF), are evaluated and compared in this study. The results demonstrate that
the RBF model exhibits higher prediction accuracy and robustness in the DFLS
application. Based on the trained RBF model, the proposed AMPC approach
constructs an associated network that directly outputs the LPV model
parameters as an adaptive, robust, and efficient prediction model. The
efficiency of the proposed approach is demonstrated through numerical
simulations of a vertical take-off thrust preparation process for the DFLS.
The simulation results indicate that the proposed AMPC method can effectively
control the DFLS thrust with a relative error below 3.5%.

###### keywords:

adaptive model predictive control , radial basis functions , linear parameter
varying model , ducted fan lift system , two-stroke piston engine control

[a] organization = School of Aerospace Engineering, Engineering Campus,
Universiti Sains Malaysia, city = Nibong Tebal, postcode = 14300, country =
Malaysia,

##  1 Introduction

The ducted fan lift system (DFLS) is widely used in current vertical takeoff
and landing (VTOL) aircraft, many of which use fuel engines as the ducted fan
drive devices[1, 2, 3, 4]. Compared to an electrically powered DFLS, a fuel-
engine-powered DFLS is more likely to meet the requirements of high power and
high energy density, thereby delivering superior flight performance. Two-
stroke aviation piston engines, which serve as the power unit of the DFLS,
have rapid, highly nonlinear dynamics with state and input constraints [5, 6].
In addition, the complex geometry of the ducted fan of the DFLS makes it
difficult to analyze its aerodynamic properties [7, 8]. These characteristics
make the engine-driven DFLS a multivariable system with tightly coupled
nonlinear dynamics, posing modeling and control challenges.

Many current spark ignition (SI) engines employ feed-forward control based on
a state observer and a proportional-integral (PI) type feedback control[9,
10]. Typically, look-up tables are used to implement the PI controller, which
necessitates a laborious process of calibration and tuning. When the state of
the engine changes rapidly, control accuracy tends to decrease[10]. To address
these challenges, researchers have developed advanced control strategies [11,
12, 13] that can be applied to SI engines, enabling more precise and energy-
efficient control. For instance, a global optimal control method based on
$H\infty$ theory was proposed for systematic control of air-fuel ratio (AFR)
with high robustness and quick responses [12]. Experimental results applied to
a four-cylinder multi-port injection (MPI) engine indicate that the AFR
control error can be limited to within 3% across a broad spectrum of operating
conditions. However, the sixth-order controller is unsuitable for real-time
computation, restricting its practical applications. Another approach applied
the dynamical sliding mode control (SMC) with a radial basis function (RBF)
neural network model to a piston engine [13]. This work demonstrated that the
SMC algorithm is robust, fast, and insensitive to parameter changes and
external disturbances within the context of nonlinear system control problems.
However, it’s important to acknowledge that as the system state approaches the
sliding mode surface, achieving precise sliding along the surface to reach the
equilibrium point becomes challenging, resulting in non-convergence.

In contrast, Model predictive control (MPC) is widely recognized as an
advanced control method in practical control engineering[14, 15]. It can
effectively address multivariable constrained optimal control problems and
offers the advantages of simplicity, straightforward design, high stability,
robustness, and adaptability[16, 17, 18]. MPC has proven to be effective for
fuel engine systems containing multiple variables, nonlinear dynamics, and
time delays[19, 20]. However, the dynamics of the DFLS can be highly
nonlinear, particularly concerning the throttle position (TPS) and the
injection fuel mass flow. It may even exhibit open-loop instability during the
transition between operation points. Thus, the performance of a linear MPC
designed for a specific operating condition will degrade near another
operating point. Recent studies[21, 22, 23] explored gain-scheduled MPC to
cope with nonlinear systems when the plant models have different orders or
time delays. This approach involves incorporating multiple predictive
controllers into a gain-scheduled MPC for various operating points, switchable
based on a predetermined scheduling policy. However, gain-scheduled MPC faces
limitations in ensuring control accuracy during transitions for nonlinear
systems and demands substantial computing resources [21]. In the realm of SI
engine control [24, 19], nonlinear MPC has gained traction due to its
potential for more precise control performance. Yet, achieving effective and
stable control performance necessitates meticulous considerations, such as
computational complexity, convergence challenges, and parameter selection,
particularly in real-time or large-scale applications.

In piston engine control, adaptive Model Predictive Control (AMPC) has emerged
as an alternative, utilizing an online-identified linear model through
successive linearization or online model estimation. This approach offers
enhanced efficiency compared to nonlinear MPC[25, 13, 26, 20]. The successive
linearization method employs a set of nonlinear ordinary differential and
algebraic equations to build the plant model and derives the linear time-
invariant (LTI) approximation at the current operating condition to update the
model parameters. However, when dealing with highly nonlinear systems,
successive linearization techniques typically require a large number of
iterations and computational resources. Online AMPC, on the other hand, has
good control accuracy and robustness, but it still comes at the cost of
significantly increased control optimization time and memory requirements due
to extensive online operations. During the data sample identification process,
online AMPC is susceptible to model errors stemming from issues such as
delays, noise, and insufficient excitation.

The objective of this paper is to develop a prediction model for AMPC that
mitigates the need for complex and time-consuming online operations. Our
proposed approach involves the direct construction of an associated linear
parameter varying (LPV) model, derived from a nonlinear global model. To
achieve this, the global model will be approximated using neural networks,
chosen for their high approximation capability and successful applications in
modeling for predictive engine control[26, 20, 27]. This paper investigated
and compared the performance of three different network models: the multi-
layer perceptron (MLP) [28], radial basis function (RBF) [29], and the Elman
network [30]. Based on the experimental results, the RBF network emerges as
the most suitable choice for representing the SI engine in the context of DFLS
AMPC. By deriving an associated linear parameter varying (LPV) model from the
global network model, we depart from the conventional practice of employing
the nonlinear model network as a prediction model or identifying a linear
model online. This shift reduces runtime computational demands and memory
usage while simultaneously bolstering DFLS control robustness. In this paper,
the network model undergoes training using the input and output data generated
by a mean value engine model (MVEM)[31], and the ducted fan model of the DFLS
will be derived using the method of theoretical design and evaluation[32].
This study focuses on the basic thrust control of the DFLS, while the thrust
fine-tuning and direction are controlled by the exit louvers configured at the
outlets of the duct[33]. Specifically, the control strategy of the DFLS
optimizes and constrains the working state of the engine by establishing a
desired thrust baseline engaged by the engine output power and a desired AFR.
In addition, based on the multistep-ahead prediction of both the aerodynamic
thrust and the AFR, the optimal control to track and maintain the desired
value is obtained as the engine state changes.

This paper presents the following main contributions: 1) This study
investigates and compares the MLP, Elman, and RBF network models, and finds
that the RBF model is more accurate in its predictions and more robust in the
DFLS application. 2) The study proposes an innovative AMPC method employing a
novel associated LPV model directly derived from an off-line trained RBF
global model, to efficiently update the prediction model state. 3) The
proposed RBF model-based AMPC is applied to the nonlinear control of a DFLS,
demonstrating the practicality of the method.

The remainder of the paper is structured as follows: Section 2 introduces the
DFLS dynamic model consisting of an MVEM and a theoretical ducted fan model.
In Section 3, three types of engine neural network models are established,
trained, and compared. Section 4 proposes the AMPC controller with an LPV
model generated from the model network trained in Section 3. Section 5
concludes the paper and outlines future research directions.

##  2 The ducted fan lift system

The DFLS dynamic model is comprised of the SI engine and ducted fan modules.
As the propulsion system, the SI engine utilizes the transmission mechanism to
finally rotate the fan shaft and generate lift force. Figure1 depicts the
control test bench scheme of the DFLS demonstrator, illustrating the assembly
relationship between the components.

Figure 1: The DFLS control test bench scheme (1: ducted fan, 2: installation
rack for the ducted fan, 3: torque and rotational speed measuring unit, 4:
belt pulley, 5: forcemeter, 6: sliding rail, 7: pulley belt, 8: SI engine, 9:
belt pulley, 10: test bench base).

###  2.1 The engine’s dynamics

For the development of engine controllers and their validation with
simulations, mathematical modeling of engine dynamics is essential[34]. The
MVEM is a widely used mathematical engine model that has achieved many
successes in real-time simulation and control of automobile and ship
engines[35, 36, 37]. The general MVEM uses empirical equations to construct
engine sub-models, which reduces engine modeling time and computational cost
by a significant margin. The implementation of the MVEM combines quasi-static
and volumetric models, dividing the two-stroke engine into four independent
volumetric control units[31]: the dynamics of intake manifold, the crankcase
and cylinder module, the dynamics of fuel injection, and the dynamics of the
crankshaft.

The equation for the crankshaft speed state in the general MVEM is[31]

$\dot{\omega}=\frac{{H}_{u}\eta_{i}(1-k_{f})\dot{m}_{f}(t-\tau_{d})}{I\omega}-\frac{{P}_{f}+{P}_{b}}{I\omega},$
(1)

where $\dot{\omega}$ is the angular acceleration of the crankshaft, $H_{u}$ is
the lower heating value of the fuel, $\eta_{i}$ is the thermal efficiency, $t$
is the time, $\tau_{d}$ is the injection-torque time delay, and $k_{f}$ is the
proportionality coefficient constant for fuel loss resulting from short-
circuiting and overflow losses during the scavenging process. $I$ is the
inertia of the engine, $P_{f}$ is the friction loss, and $P_{b}$ is the load
power of the engine. The torque of the engine can be expressed as[31]

${Q}_{eng}=\frac{{H}_{u}\eta_{i}(1-k_{f})\dot{m}_{f}(t-\tau_{d})}{\omega}-\frac{{P}_{f}}{\omega}.$
(2)

The intake airflow is mixed with fuel during the engine’s intake process, and
the normalized air-fuel ratio is

$\lambda={\frac{\dot{m}_{as}}{\dot{m}_{f}L_{th}}},$ (3)

where $L_{th}$ is the stoichiometric air-fuel ratio.

###  2.2 Ducted fan dynamics

Due to the complexity of the geometry, it is difficult to analyze the
aerodynamic properties of ducted fans. This paper adopts a modeling method for
the DFLS based on the blade element theory (BET) and the momentum theory[32].
According to the BET, propellers are made up of minuscule elements in the
shape of airfoils along the radius of each blade[38]. The resulting velocity
for each element can be decomposed into rotational and translational
components. The resulting aerodynamic force can be decomposed into the drag
and lift. When decomposed in the plane of rotation, it produces thrust and the
torque-producing force.

The total thrust $T_{UDF}$ can be expressed as[32]

$T_{UDF}={\frac{1}{2}}{\rho}{V_{trans}}^{2}B\int_{0}^{R}{T_{c}\cdot{dr}},$ (4)

where $\rho$ is the air density, $V_{trans}$ is the translational speed
component, $B$ is the blade number corrective factor, $R$ is the blade radius,
$T_{c}$ is the thrust coefficient and $dr$ represents the infinitesimal
airfoils along the blade radius. Similarly, the total torque $Q_{UDF}$ can be
expressed as:

$Q_{UDF}={\frac{1}{2}}{\rho}{V_{trans}}^{2}B\int_{0}^{R}{Q_{c}\cdot{dr}},$ (5)

where $Q_{c}$ is the torque coefficient. The rate of energy supplied by the
engine matches its power $P_{UDF}$:

$P_{UDF}=2\pi{n}{Q_{UDF}}.$ (6)

Using the momentum theorem, the static thrust ratio between the ducted fan and
the unducted fan with the same absorbed power[39] can be calculated as:

${\frac{T_{DF}}{T_{UDF}}}=1.26\Bigg{(}{\frac{S_{3}}{S_{2}}}\Bigg{)}^{\frac{1}{3}},$
(7)

where $S_{2}$ is the area of the fan disc and $S_{3}$ is the area of the duct
outlet. Equations 1, 2, 4-7 can be used to calculate the thrust of the ducted
fan that corresponds to the engine output in this instance. Specifically, the
absorbed power of the ducted fan is determined by calculating the engine
output power from the MVEM output. BET is then used to calculate the thrust of
an unducted fan with the same absorption power.

Previous studies have validated the dynamic models of the DFLS[31, 32]. Both
models are explicit mathematical representations, allowing for flexible
combination and application. In this study, the established DFLS dynamic model
will be utilized as the plant for control research, and the MVEM data will be
employed to create a neural network engine model intended for the network-
based controller in the following section.

##  3 Neural network model of the engine

As stated previously, the engine network model is established using MLP,
Elman, and RBF neural networks in this study. After comparing the accuracy and
robustness of these three network models, the best one will be utilized in the
development of the AMPC controller. Figure 2 depicts the expanded engine
model, which has four inputs, fuel injection rate $\dot{m}_{fi}$, throttle
position $TPS$, engine speed $n$, and normalized AFR $\lambda$, and generates
three outputs, normalized AFR $\lambda$, engine speed $n$, and engine output
torque $Q_{eng}$ at next instance.

Figure 2: Expanded engine model structure

The engine network model developed will not be employed as an adaptive
prediction model for the MPC directly, but rather to facilitate the
development of an LPV prediction model generated from the trained network
parameters.

###  3.1 MLP network model

Based on the gradient descent algorithm, the MLP network is a supervised
learning technique[40]. The MLP network structure is depicted in Figure 3 and
the mathematical expression of the MLP network output is

Figure 3: The concise structure of the MLP network.

$\bm{h}=f_{2}(LW\cdot f_{1}(IW\cdot\bm{p}^{T}+\bm{b}_{1})+\bm{b}_{2}),$ (8)

where $\bm{p}$ represents the input vector, $IW$ represents the weight matrix
from the input layer to the hidden layer, and $LW$ represents the weight
matrix from the hidden layer to the output layer. The bias vectors for the
hidden layer and output layer are denoted by $\bm{b}_{1}$ and $\bm{b}_{2}$,
respectively.

The primary learning phases of MLP networks are mode forward transmission
process and error back propagation (BP). In the back propagation process, when
the output does not achieve the desired value, the error signal transfers back
along the original path. Adjustment of general functions based on cumulative
value and domain is expressed as[41]:

$W(k+1)=W(k)-\alpha\frac{\partial E_{k}}{\partial W(k)},$ (9)
$\bm{b}(k+1)=\bm{b}(k)-\beta\frac{\partial E_{k}}{\partial b(k)},$ (10)

where $W(k+1)$ and $W(k)$ are the weights at times $k+1$ and $k$ respectively,
and $\bm{b}(k+1)$ and $\bm{b}(k)$ are the biases at times $k+1$ and $k$
respectively. The error function is denoted by $E_{k}$, and $\alpha$ and
$\beta$ are the learning rates.

For the MLP network, different number of hidden layer nodes have been used in
training experiments and a structure with a single hidden layer containing 26
nodes is chosen as shown in Figure 2, which gives the minimum prediction
error. The activation function of the hidden layer $f_{1}$ is hyperbolic
tangent,and the activation function of the output layer $f_{2}$ is linear
transfers. The learning rates $\alpha$ and $\beta$ are set as 0.1 in this MLP
network to ensure the convergence of network training. In this paper, the
training accuracy target is set to 0.0001.

###  3.2 Elman network model

Elman network is a well-known partial recurrent network, which lies between a
traditional feed-forward perception network and a pure recurrent network[30].
As shown in Figure 4, a back-forward loop is included in an Elman network that
is sensitive to input data history[42].

Figure 4: The structure of the Elman network.

The Elman network is also trained using the dynamic BP algorithm, and its
hidden layer activation function $f_{1}$ utilizes the hyperbolic tangent
function, and the output layer activation function $f_{2}$ utilizes the linear
transfer function. Due to the additional backforward loop, the Elman network’s
input-output relations differ from those of the MLP network, as illustrated
below:

$\textbf{a}(k)=f_{1}(IW\cdot\textbf{p}+LW_{1}\cdot\textbf{a}(k-1)+\textbf{b}_{1}),$
(11) $\textbf{h}(k)=f_{2}(LW_{2}\cdot\textbf{a}(k)+\textbf{b}_{2}).$ (12)

The Elman network selects a structure with a single hidden layer which
contains 12 nodes, after the training experiments with different number of
hidden layer nodes. The learning rate is set to be 0.01, and the training
accuracy target is set to 0.0001. The maximum number of training epochs is set
as 1000.

###  3.3 RBF network model

RBF neural networks consist of an input layer, a single hidden layer with a
radial basis activation function, and a linear output layer[42]. Gaussian is
the most frequently used activation function:

$\phi(\textbf{{p}},\textbf{{c}})=e^{-{\frac{\|\textbf{{p}}-\textbf{{c}}\|^{2}}{s^{2}}}},$
(13)

where $\bm{c}$ is the center of the Gaussian function and $s$ is the radius,
which gives a measure of the spread of the Gaussian curve. Figure 5 depicts
the architecture of the RBF network. The distance between the input vector
$\bm{p}$ and the center vector $\bm{c}$ is denoted by $\|dist\|$, and the RBF
network output is:

$\bm{h}=LW\cdot e^{-{\frac{\|\textbf{{p}}-\textbf{{c}}\|^{2}}{s^{2}}}}.$ (14)

Figure 5: The concise structure of the RBF network.

The learning process of the RBF network can be divided into three stages: in
the first stage, based on the distribution of the input sample, the centers
and radius values of each node in the hidden layer are determined. In the
second stage, the output layer weights are calculated using least-square
methods. Thirdly, the parameters of the hidden layer and output layer are
simultaneously adjusted based on the sample signal to improve the network
precision. Using the Poggio[43] method and the r-nearest neighborhood
heuristic, the authors determined the centers $c$ and the radius $s$ of the
hidden layer nodes in the RBF network. For training the network weights, the
least-mean-square (LMS) algorithm is used. In this study, various numbers of
hidden layer nodes have been tested, and a hidden layer with 25 centers is
selected for the RBF network.

###  3.4 Training and comparison of neural network models

In the engine data collection stage, training data must be representative of
typical plant behaviors in order to evaluate the performance of various engine
models under practical driving conditions. Consequently, the sampled data
should adequately represent the state space of the system to be controlled.
The time scale of the MVEM is just sufficient to accurately describe the
changes in the engine variables that change the most rapidly, which is
advantageous for engine control applications [44]. The trained neural network
model will therefore have adequate transient and steady-state performance.

A total of 1000 data samples are generated and separated into two groups: 950
samples for training and 50 samples for validation. Before the training and
validation, all input and output data were normalized to the range of
$[-1,1]$. The Gaussian white noise (GWN) is added to the neural network
training samples and its signal-to-noise ratio (SNR) is 5 db. Figures 6-8
compare the prediction results and the corresponding proportional error (PE)
over the validation data after training with 950 samples.

a) Prediction results from network models

b) PE of prediction results

Figure 6: Comparison of predicted engine torque.

a) Prediction results from network models

b) PE of prediction results

Figure 7: Comparison of predicted engine speed.

a) Prediction results from network models

b) PE of prediction results

Figure 8: Comparison of predicted engine AFR.

The throttle position is constrained between 5% and 90%, and the fuel
injection rate is between 0.0011 and 0.0055 kg/s. The sampling interval has
been set to 0.1 s. Table 1 displays the mean absolute percentage error (MAPE)
of the prediction results from the three network models.

Table 1: The MAPE of the prediction results. | MLP | RBF | Elman
---|---|---|---
Engine torque | 0.81% | 0.49% | 1.44%
Engine speed | 1.08% | 0.20% | 1.62%
AFR | 0.95% | 1.13% | 1.92%

The comparison results reveal that the prediction of the network models
correspond well to the engine output during the model validation phase, where
the maximum MAPE is 1.92%. The RBF network model has the highest accuracy with
the MAPEs of 0.49%, 0.20% and 1.13%. This is followed by the MLP model and the
Elman model, in order of accuracy. Due to the accuracy and robustness
advantages of the RBF network model, it will be used in the development of the
AMPC in the following section.

##  4 RBF-based adaptive model predictive control

Adapting the prediction model to changing operating conditions allows the AMPC
to address highly nonlinear control problems. As previously mentioned,
updating the prediction model states with an off-line LPV model is expected to
enable efficient and robust AMPC strategies. Consequently, this section
establishes an LPV model based on the parameters of the previously selected
RBF network model, facilitating the development of AMPC for the DFLS.
Specifically, the AMPC system structure is introduced first, followed by the
presentation of the LPV model constructed through a network associated with
the RBF model. Finally, simulation studies on the RBF-based AMPC for DFLS
control are conducted.

###  4.1 Control system structure

The idea of adaptive model predictive control has been introduced in details
in the literature[20]. Figure 9 shows the block diagram of the AMPC structure.
The nonlinear DFLS model is consisting of the MVEM and the ducted fan modules
in the control simulations. At each time step, the LPV model associated with
the RBF network will generate a constant linear prediction model based on the
current state and input.

Figure 9: The block diagram of the AMPC structure.

The tracking control target is to minimize the error between the system output
$\bm{y}=[T_{DF},\lambda]^{T}$ and the reference $\bm{y}^{ref}$ over a
specified time horizon. The AMPC controller will be used to optimize two DFLS
inputs, $\bm{u}=[TPS,\dot{m}_{fi}]^{T}$, to achieve precise dynamic tracking
of the desired DFLS thrust ${T_{DF}}^{ref}$ and the engine AFR
$\lambda^{ref}$. A cost function $Z$ is therefore defined as a quadratic
function in terms of tracking errors and input increments:

$\displaystyle Z(k)$
$\displaystyle=\epsilon\sum_{j=k+N_{1}}^{k+N_{2}}\\{[\lambda^{ref}(j)-\hat{\lambda}(j)]^{2}+[{T_{DF}}^{ref}(j)-\hat{T}_{DF}(j)]^{2}\\}$
(15)
$\displaystyle+\xi\sum_{j=k}^{k+N_{c}}\\{[\dot{m}_{fi}(j)-\dot{m}_{fi}(j-1)]^{2}+[TPS(j)-TPS(j-1)]^{2}\\},$

where $k$ represents the current time instance, $\hat{\lambda}$ is the
predicted AFR, $\hat{T}_{DF}$ is the predicted thrust of the DFLS, and
$\epsilon$ and $\xi$ are control weighting factors that penalize the excessive
modification of control inputs. Future horizon $N_{1}$ and $N_{2}$, also known
as the prediction horizon, specifies the number of upcoming samples. Likewise,
the control horizon $N_{c}$ specifies the number of samples from which optimal
inputs are calculated. At each sampling instant, the optimization yields a
sequence of input signals, but only the first input is applied to the
plant[45, 46] resulting in a receding horizon approach[47].

###  4.2 The LPV model

The estimation of the state and the update of the prediction model are
indispensable to the effectiveness, precision, and robustness of an AMPC[26].
The AMPC has currently adopted three categories of model updating strategies:
online parameter estimation, successive linearization, and linear parameter
varying (LPV) modeling. The online parameter estimation method estimates model
parameters from real-time plant measurements and updates the prediction model.
However, it can be computationally intensive, and suitable for control
applications with longer control intervals and sufficient resources. The
successive linearization method builds the plant model using nonlinear
equations and derives a linear approximation to update the model parameters.
But dealing with highly nonlinear systems using successive linearization often
requires many iterations and computations. An LPV system is a linear state-
space model whose dynamics vary as a function of certain time-varying
parameters. Using an off-line LPV model has the advantage of supporting rapid
batch linearization to obtain a variety of plant models at the desired
operating points for model updating. This is advantageous for the computation
efficiency and stability of the DFLS AMPC. However, LPV models of nonlinear
mechanical systems are typically constructed based on the system’s dynamic
model, which is challenging for SI engines.

This paper proposes to directly utilize the trained RBF network to construct
the LPV model for the AMPC prediction step. Specifically, an associated
network can be derived from the RBF model off-line[48] and outputs directly
the Jacobian matrix of the RBF model, which contains the LPV model parameters.
The LPV model of the DFLS can be expressed in the discrete linear form with
time-varying parameters:

$\Delta\bm{x}(t+1)=A(t)\Delta\bm{x}(t)+B(t)\Delta\bm{u}(t),$ (16)
$\Delta\bm{y}(t)=C(t)\Delta\bm{x}(t)+D(t)\Delta\bm{u}(t),$ (17)

where the states, inputs, and outputs of the DFLS are

$\bm{x}(t)=[Q_{eng}(t),n(t),\lambda(t)],$ (18)
$\bm{u}(t)=[TPS(t),\dot{m}_{fi}(t)],and$ (19)
$\bm{y}(t)=[T_{DF}(t),\lambda(t)].$ (20)

Because the current state $Q_{eng}(t)$ has no effect on the subsequent system
states, the first column of the system matrix $A(t)$ contains only zeros. And
the linearized system’s matrices can be written as follows:

$A(t)=\frac{\partial\bm{x}(t+1)}{\partial\bm{x}(t)}=\begin{bmatrix}0&\frac{\partial
Q_{eng}(t+1)}{\partial n(t)}&\frac{\partial
Q_{eng}(t+1)}{\partial\lambda(t)}\\\ 0&\frac{\partial n(t+1)}{\partial
n(t)}&\frac{\partial n(t+1)}{\partial\lambda(t)}\\\
0&\frac{\partial\lambda(t+1)}{\partial
n(t)}&\frac{\partial\lambda(t+1)}{\partial\lambda(t)}\end{bmatrix},$ (21)

$B(t)=\frac{\partial\bm{x}(t+1)}{\partial\bm{u}(t)}=\begin{bmatrix}\frac{\partial
Q_{eng}(t+1)}{\partial TPS(t)}&\frac{\partial
Q_{eng}(t+1)}{\partial\dot{m}_{fi}(t)}\\\ \frac{\partial n(t+1)}{\partial
TPS(t)}&\frac{\partial n(t+1)}{\partial\dot{m}_{fi}(t)}\\\
\frac{\partial\lambda(t+1)}{\partial
TPS(t)}&\frac{\partial\lambda(t+1)}{\partial\dot{m}_{fi}(t)}\end{bmatrix},$
(22)

$C(t)=\frac{\partial\bm{y}(t+1)}{\partial\bm{x}(t+1)}=\begin{bmatrix}\frac{\partial
T_{DF}(t+1)}{\partial Q_{eng}(t+1)}&\frac{\partial T_{DF}(t+1)}{\partial
n(t+1)}&0\\\ 0&0&1\end{bmatrix},$ (23)

$D(t)=\frac{\partial\bm{y}(t+1)}{\partial\bm{u}(t+1)}=\begin{bmatrix}0&0\\\
0&0\end{bmatrix}.$ (24)

The output matrix $C$ can be calculated explicitly using equations 4-7 from
the ducted fan model. And the system matrix $A$ and the control effectiveness
matrix $B$ are associated with the SI engine, the elements of which can be
approximated by the partial derivatives of the trained RBF network inputs
$\textbf{{h}}=[TPS(t),\dot{m}_{fi}(t),n(t),\lambda(t)]$ with respect to the
network outputs $\textbf{{p}}=[Q_{eng}(t+1),n(t+1),\lambda(t+1)]$.

To determine the Jacobian matrix of the RBF model $\mathcal{J}$, we can
rewrite the $j$th Gaussian radial function in the RBF network from equation 13
as below:

$\phi_{j}(p)=e^{-({\|\bm{p}-\bm{c}_{j}\|}/{{s_{j}}})^{2}}.$ (25)

where $c_{j}$ is the $j$th center point and $s_{j}$ is its radius. The
associated network is linear to the parameters $LW$, which are weights
connecting the output of each radial function to the output of the network,
when $c_{j}$ and $s_{j}$ are fixed. Derivatives can therefore be explicitly
calculated as a function of the network input p and network weights $LW$:

$\mathcal{J}(\textbf{{p}},LW)=LW^{T}\Phi(\textbf{{p}})=LW^{T}\begin{bmatrix}\phi_{1}(\textbf{{p}})\frac{-2}{{s_{1}}^{2}}(\textbf{{p}}-\textbf{{c}}_{1})\\\
\vdots\\\
\phi_{J}(\textbf{{p}})\frac{-2}{{s_{J}}^{2}}(\textbf{{p}}-\textbf{{c}}_{J})\end{bmatrix}.$
(26)

And this function can be constructed as another network associated with the
RBF model network as depicted in Figure 10, with $J$ activation functions,
each of which is the partial derivative of the $j$th Gaussian radial function
output with respect to the critic input:

Figure 10: The system model network with RBF activation functions and its
associated network. These two networks share the same set of centers, denoted
by $c_{j}$, and the to-be-determined parameters $LW$.[48]

$\displaystyle\Phi_{j}(\bm{p})$
$\displaystyle=\frac{\phi_{j}(\bm{p})}{\partial\bm{p}}=e^{-({\|\bm{p}-\bm{c}_{j}\|}/{{s_{j}}})^{2}}\frac{-2}{{s_{j}}^{2}}(\bm{p}-\bm{c}_{j})$
(27)
$\displaystyle=\phi_{j}(\bm{p})\frac{-2}{{s_{j}}^{2}}(\bm{p}-\bm{c}_{j}),$

where $\Phi_{j}$ represents the $j$th row vector of the derivatives
$\Phi(\bm{p})=\partial\phi(\bm{p})/\partial\bm{p}$.

This associated network model outputs directly the Jacobian matrix of the RBF
model $\mathcal{J}$, whose elements are the partial derivatives of the model
outputs with respect to the model inputs, and which can be used to construct
the LPV model matrices $A$ and $B$ in equations 21 and 22. Note that this
concise mathematical relationship is derived from the linear-in-parameter
property of the system model [48], which will effectively enhance the AMPC
prediction model’s update speed. Moreover, employing a linear-in-parameter
system model can aid in avoiding the local minimum trap. The system model
network is intended to be nonlinear in inputs but linear in parameters, and
the associated network can be derived explicitly and precisely from the system
model network. This characteristic is indispensable for the computational
efficiency of the RBF-based AMPC.

###  4.3 AMPC simulation of the DFLS

To validate the proposed AMPC, the control simulation of the DFLS during
takeoff thrust preparation is implemented. During the vertical take-off power
preparation procedure, the anticipated thrust of the DFLS ${T_{DF}}^{ref}$
gradually increases from the idling rating (10 $kgf$) to the hovering thrust
(80 $kgf$). During this procedure, the engine gradually increases its throttle
to a relatively stable state with a low expected normalized AFR (0.82),
ensuring adequate power output, before adjusting the expected normalized AFR
(1.0) to an efficient mode. In the simulation study, the throttle position is
constrained between 5% and 90%, and the fuel injection rate is between 0.0011
and 0.0055 kg/s. The thrust of the DFLS is constrained between 0 and 150 kgf,
and the normalized AFR is between 0.68 and 1.26.

In order to compare the AMPC against a traditional MPC, a linear MPC is
designed for the same control process for the DFLS. The linear MPC and the
AMPC are designed with the same cost function, as illustrated in equation 15,
and same parameters. The nonlinear optimization parameters were set as
follows: $N_{1}=1$, $N_{2}=8$, $N_{c}=3$, $\epsilon=0.8$, and $\xi=0.5$. In
the control simulation, the normalized AFR measurement and the DFLS thrust
measurement were subjected to Gaussian noise $\mathcal{N}$(0, 0.005). Figures
11-12 depict the simulation tracking results for the DFLS’s desired thrust and
AFR. Before time step $100$, the desired thrust is continuously increasing,
the tracing control results of the DFLS thrust and the engine AFR with the
linear MPC are divergent.

Figure 11: The DFLS thrust tracking performance with a linear MPC. Figure 12:
The AFR tracking performance with a linear MPC.

After the expected thrust has been stabilized, the linear MPC can only
precisely track the expected thrust, while the engine AFR control result is
inconsistent with the desired value. As expected, the linear MPC’s tracking
and regulating performances are unacceptable.

The AMPC parameters are manually adjusted to determine the appropriate values
by comparing the control performance. Figure 13 depicts the tracking control
results for the desired thrust. Figure 14 depicts the corresponding simulation
results of the tracking curve for the desired normalized AFR. Figure 15
illustrates the AMPC-optimized $TPS$ and $\dot{m}_{fi}$ engine results,
respectively. Figure 16 depicts the corresponding simulation results for the
engine speed and output torque.

a) Tracking performance of the desired thrust of the DFLS

b) Relative error of the thrust control

Figure 13: The DFLS thrust tracking performance with the proposed AMPC.

a) Tracking performance of the desired AFR

b) Relative error of the AFR control

Figure 14: The AFR tracking performance with the proposed AMPC.

a) Optimized throttle position by the AMPC

b) Optimized fuel injection mass flow by the AMPC

Figure 15: The control inputs of the engine (the TPS and the fuel injection
mass flow $\dot{m}_{fi}$) optimized by the proposed AMPC.

a) The engine crankshaft speed

b) The engine crankshaft torque

Figure 16: The engine crankshaft speed $n$ and the engine crankshaft torque
$Q_{eng}$ during the control simulation.

The simulation results indicate that the developed AMPC can accurately control
the thrust and AFR during the vertical take-off power preparation process. As
illustrated in Figures 13,14, the relative error ranges for DFLS thrust and SI
engine AFR tracking control are -3.5% to 1.0% and -2.1% to 2.2%, respectively.
The aforementioned outcomes validated the effectiveness and robustness of the
proposed RBF-based AMPC.

##  5 Conclusions

This paper presents a novel adaptive model predictive control (AMPC) approach
for controlling the thrust of an engine-driven ducted fan lift system (DFLS).
The proposed method is based on an off-line linear parameter varying (LPV)
model, whose parameter updating laws is derived from a radial basis function
(RBF) network. The use of the RBF network is motivated by its superior
prediction accuracy and robustness compared to other network models, such as
multi-Layer perceptron and Elman. This LPV model enables real-time updates of
the AMPC prediction model across the full operating envelope, making it an
effective solution for handling the highly nonlinear dynamics of DFLS. The
proposed AMPC receives and updates its LPV parameters from an associated
network without massive online operations, which enhances the control
effectiveness and avoids model errors resulting from issues such as delays,
noise, and insufficient excitation. This concise mathematical relationship is
derived from the linear-in-parameter property of the system model, avoiding
local minimum traps. The DFLS AMPC was validated in numerical simulations,
demonstrating its ability to achieve precise control of thrust and air-fuel
ratio during the vertical take-off preparation process. The control strategy
is designed to track the desired thrust by controlling the engine output power
while ensuring reliability and efficiency through synchronized air-fuel ratio
control.

The RBF model-based AMPC approach proposed in this paper is efficient and
practical. The validation results show its potential for wider application to
other nonlinear industrial control problems. However, the LPV model obtained
off-line may lack sufficient robustness to handle stochastic and uncertain
plant model changes. Further research is needed to explore the combination of
the proposed method with online AMPC to preserve control efficiency and
enhance its adaptability. In conclusion, this study presents a promising new
method for controlling the thrust of engine-driven ducted fan lift systems and
opens up avenues for further improvement and refinement.

## Acknowledgement

The corresponding author would like to thank Malaysian Ministry of Higher
Education (MOHE) for providing the Fundamental Research Grant Scheme (FRGS):
FRGS/1/2020/TK0/USM/03/11.

## References

  * [1] G. Warwick.  Boeing to study ducted-fan vtol with israel’s tactical robotics.  Aerospace Daily & Defense Report, (15):271, 2020.
  * [2] M. Speck, Jrg Buchholz, and M. Sellier.  A mathematical model of a twin ducted fan vtol jetpack.  Proceedings of the Institution of Mechanical Engineers Part G Journal of Aerospace Engineering, 1(1):1–47, 2015.
  * [3] P. Host.  Martin uav, usmc working on v-bat hybrid-electric, heavy fuel propulsion.  Jane s international defense review: IHS Jane s international defense review, (Oct.):54, 2021.
  * [4] Frank and Colucci.  Darpa’s aerial reconfigurable embedded system(ares) leverages ducted fan technology to resupply small combat units in tight spots.  Vertiflite, 62(5):28–31, 2016.
  * [5] A. Balluchi, L. Benvenuti, D. Benedetto, C. Pinello, and A. L. Sangiovanni-Vincentelli.  Automotive engine control and hybrid systems: challenges and opportunities.  Proc IEEE, 88(7):888–912, 2000.
  * [6] Yixuan, WANG, Yan, SHI, Maolin, CAI, Weiqing, XU, Qihui, and YU.  Efficiency optimized fuel supply strategy of aircraft engine based on air-fuel ratio control.  Chinese Journal of Aeronautics, 32(2):489–498, 2019.
  * [7] T. Ai, B. Xu, C. Xiang, W. Fan, and Y. Zhang.  Aerodynamic analysis and control for a novel coaxial ducted fan aerial robot in ground effect:.  International Journal of Advanced Robotic Systems, 17(4):1–26, 2020.
  * [8] N.-J Fan, J. Xu, and J.-W Yang.  Numerical study on aerodynamic characteristics of ducted fan aircraft.  Journal of Beijing Insititute of Technology, 19:19–24, 03 2010.
  * [9] Munan Hong, Minggao Ouyang, Tielong Shen, and Jianqiu Li.  Model-based pi feedback control of engine torque.  In 2010 8th IEEE International Conference on Control and Automation. [v.1], pages 12–15, Xiamen, China, 2010. Se Key Lborory of Auooive Sfey nd Energy, Tsinghu Universiy.
  * [10] M. Won, S. B. Choi, and J. K. Hedrick.  Air-to-fuel ratio control of spark ignition engines using gaussian network sliding control.  IEEE Transactions on Control Systems Technology, 6(5):678–687, 2002.
  * [11] Ye Zhou, Hann Woei Ho, and Qiping Chu.  Extended incremental nonlinear dynamic inversion for optical flow control of micro air vehicles.  Aerospace Science and Technology, 116:106889, 2021.
  * [12] Christian Vigild, Karsten Andersen, Elbert Hendricks, and Michael Struwe.  Towards robust h-infinity control of an si engine’s air/fuel ratio.  SAE 1999 Transactions - Journal of Engines, pages 1–16, 03 1999.
  * [13] S. Wang and D. L. Yu.  Adaptive rbf network for parameter estimation and stable air-fuel ratio control.  Neural Networks, 21(1):102–122, 2008.
  * [14] Y. U. Bing, L. I. Zhouyang, K. E. Hongwei, and T. Zhang.  Wide-range model predictive control for aero-engine transient state.  Chinese Journal of Aeronautics, 35(7):246–260, 2022.
  * [15] Guillaume J. J. Ducard Leonard Bauersfeld, Lukas Spannagl and Christopher H. Onder.  Mpc flight control for a tilt-rotor vtol aircraft.  IEEE Transactions on Aerospace and Electronic Systems, 57(4):2395 – 2409, 2021.
  * [16] Eduardo Camacho and Carlos Bordons.  Model predictive control.  Springer London, 01 2013.
  * [17] Kaian Chen, Kaixiang Zhang, Zhaojian Li, Yan Wang, Kai Wu, and V Kalabi, Uros.  Stochastic model predictive control for quasi-linear parameter varying systems: Case study on automotive engine control.  Journal of Dynamic Systems, Measurement, and Control, (144-6):1–9, 2022.
  * [18] Gokul Siva Sankar, Rohan Shekhar, C. Manzie, Takeshi Sano, and Hayato Nakada.  Fast calibration of a robust model predictive controller for diesel engine airpath.  IEEE Transactions on Control Systems Technology, 28:1–16, 04 2018.
  * [19] Na Li, Wei Gao, and Yunfeng Hu.  Nonlinear model predictive controller design for air system control of variable geometry turbocharged diesel engine.  IFAC-PapersOnLine, 51(31):714–719, 2018.  5th IFAC Conference on Engine and Powertrain Control, Simulation and Modeling E-COSM 2018.
  * [20] Samir Saraswati and Satish Chand.  Online linearization-based neural predictive control of air–fuel ratio in si engines with pid feedback correction scheme.  Neural Computing & Applications, 19(6):919–933, 2010.
  * [21] J. Gao and H. M. Budman.  Design of robust gain-scheduled mpc controllers for nonlinear processes.  IFAC Proceedings Volumes, 39(2):323–328, 2006.
  * [22] A. Ilka and V. Vesely.  Gain-scheduled mpc design for nonlinear systems with input constraints.  Ifac Papersonline, 48(11):912–917, 2015.
  * [23] L.P. Carvalho, J.M. Palma, T.E. Rosa, B. Jayawardhana, and O.L.V. Costa.  Gain-scheduled controller for fault accommodation in linear parameter varying systems with imprecise measurements.  IFAC-PapersOnLine, 54(8):57–63, 2021.  4th IFAC Workshop on Linear Parameter Varying Systems LPVS 2021.
  * [24] Kazuyoshi, Odachi, Hiroyuki, Oyama, Masaki, and Yamakita.  Nonlinear model predictive control of si engine using discrete time identification model.  In The Society of Instrument and Control Engineers Annual Conference 2015, volume 4, pages 51–56, 2015.
  * [25] Yu, Feng, Xiaohong, Jiao, Zhijing, and Wang.  Adaptive and predictive control of fuel injection with time delay consideration for air–fuel ratio regulation of gasoline direct injection engines.  Proceedings of the Institution of Mechanical Engineers Part D Journal of Automobile Engineering, pages 1–12, 2017.
  * [26] Y. J. Zhai and D. L. Yu.  Neural network model-based automotive engine air/fuel ratio control and robustness evaluation.  Engineering Applications of Artificial Intelligence, 22(2):171–180, 2009.
  * [27] J. Lu, Z. Xu, and L. Liu.  Compression ratio control of an opposed-piston free-piston engine generator based on artificial neural networks.  IEEE Access, PP(99):107865–107875, 2020.
  * [28] Q. Zhao and T. Higuchi.  Efficient learning of nn-mlp based on individual evolutionary algorithm.  Neurocomputing, 13(2-4):201–215, 1996.
  * [29] Fma Acosta.  Radial basis function and related models: An overview.  Signal Processing, 45(1):37–58, 1995.
  * [30] A. Kalinli and S. Sagiroglu.  Elman network with embedded memory for system identification.  Journal of Information Science and Engineering, 22(6):1555–1568, 2006.
  * [31] Hanjie Jiang, Ye Zhou, Hann Woei Ho, and Elmi Abu Bakar.  Modeling of two-stroke aviation piston engines for control applications.  Advances in Mechanical Engineering, 15(2):1–15, 2023.
  * [32] Hanjie Jiang, Ye Zhou, and Hann Woei Ho.  Aerodynamic design and evaluation of a ducted fan lift system for vertical takeoff and landing flying cars.  Proceedings of the Institution of Mechanical Engineers, Part A: Journal of Power and Energy, 237(1):115–125, 2023.
  * [33] J. A. Yuska and J. H. Diedrich.  Fan and wing force data from wind tunnel investigation of a 0.38 meter (15 inch) diameter vtol model lift fan installed in a two dimensional wing.  NASA TECHNICAL NOTE, pages 1–101, 1972.
  * [34] J. N. Dias, G. O. Violato, and C. A. Martins.  Dynamic model of a two-stroke glow engine from experimental data.  Proceedings of the Institution of Mechanical Engineers Part G Journal of Aerospace Engineering, 226(12):1502–1512, 2012.
  * [35] R. Rajamani.  Mean value modeling of si and diesel engines.  Vehicle Dynamics and Control, pages pp 241–265, 2012.
  * [36] G. Theotokatos, C. Guan, H. Chen, and I. Lazakis.  Development of an extended mean value engine model for predicting the marine two-stroke engine operation at varying settings.  Energy, page S0360544217318431, 2018.
  * [37] Kevin, KoosupYum, BhushanTaskar, EilifPedersen, and SverreSteen.  Simulation of a two-stroke diesel engine for propulsion in waves.  International Journal of Naval Architecture and Ocean Engineering, 9(4):351–372, 2017.
  * [38] M. Eriksson and B. Wedell.  Performance estimation of a ducted fan uav.  Linköping University, Department of Electrical Engineering, 2006.
  * [39] Jr Platt, Robert J.  Static tests of a shrouded and an unshrouded propeller.  Technical Report Archive & Image Library, (L7H25):1–40, 02 1948.
  * [40] Tadeusz Wilusz.  Neural networks — a comprehensive foundation.  Neurocomputing, 8(3):359–360, 1995.
  * [41] F. Jia, H. Li, R. Lai, H. Yao, and Z. Bao.  Evaluation of aerospace production site management and control capability based on f-ahp and bp neural network.  IOP Conference Series: Materials Science and Engineering, 1043(4):042027 (12pp), 2021.
  * [42] X. Tong, W. Zhu, and H. Yu.  A research using hybrid rbf/elman neural networks for intrusion detection system secure model.  Computer Physics Communications, 180(10):1795–1801, 2009.
  * [43] T. Poggio.  A theory of networks for approximation and learning.  A.i.memo, 1140(9):1481 – 1497, 1989.
  * [44] E. Hendricks and Spencer C. Sorenson.  Mean value modelling of spark ignition engines.  SAE transactions, pages 1359–1373, 1990.
  * [45] Alessandro Alessio and Alberto Bemporad.  A Survey on Explicit Model Predictive Control, pages 345–369.  Springer Berlin Heidelberg, Berlin, Heidelberg, 2009.
  * [46] Jay, H., and Lee.  Model predictive control: Review of the three decades of development.  International Journal of Control Automation & Systems, 9:415–424, 2011.
  * [47] S. W. Wang, D. L. Yu, J. B. Gomm, G. F. Page, and S. S. Douglas.  Adaptive neural network model based predictive control for air-fuel ratio of si engines.  Engineering Applications of Artificial Intelligence: The International Journal of Intelligent Real-Time Automation, (2):189–200, 2006.
  * [48] Ye Zhou.  Efficient online globalized dual heuristic programming with an associated dual network.  IEEE Transactions on Neural Networks and Learning Systems, pages 1–12, 2022.


REMIRO-AZÓCAR et al

*Antonio Remiro-Azócar, Medical Affairs Statistics, Bayer plc, Reading, United Kingdom. . Tel: (+44) 7818 746 781

Antonio Remiro-Azócar, Medical Affairs Statistics, Bayer plc, 400 South Oak
Way, Reading, RG2 6AD, United Kingdom

# Methodological considerations for novel approaches to covariate-adjusted
indirect treatment comparisons

Antonio Remiro-Azócar  Anna Heath  Gianluca Baio  Medical Affairs Statistics,
Bayer plc, Reading, United Kingdom  Department of Statistical Science,
University College London, London, United Kingdom  Child Health Evaluative
Sciences, The Hospital for Sick Children, Toronto, Canada  Dalla Lana School
of Public Health, University of Toronto, Toronto, Canada  antonio.remiro-
<EMAIL_ADDRESS>
(26 April 2016; 6 June 2016; 6 June 2016)

###### Abstract

We examine four important considerations in the development of covariate
adjustment methodologies for indirect treatment comparisons. Firstly, we
consider potential advantages of weighting versus outcome modeling, placing
focus on bias-robustness. Secondly, we outline why model-based extrapolation
may be required and useful, in the specific context of indirect treatment
comparisons with limited overlap. Thirdly, we describe challenges for
covariate adjustment based on data-adaptive outcome modeling. Finally, we
offer further perspectives on the promise of doubly-robust covariate
adjustment frameworks.

###### keywords:

Health technology assessment, indirect treatment comparison, matching-adjusted
indirect comparison, parametric G-computation, covariate adjustment, model
misspecification

††articletype: Discussion

Covariate-adjusted indirect comparisons are increasingly used in health
technology assessment (HTA) to compare treatments that have not been trialed
against each other.1, 2 Matching-adjusted indirect comparison (MAIC)3 is the
most popular methodology.4 MAIC weights the subjects in one of the studies by
the odds of being assigned to another study, conditional on a set of baseline
covariates, without explicitly requiring an outcome model. The estimated
weights enforce balance in selected covariate moments between studies.

Recently, alternatives based on estimating the conditional outcome expectation
given treatment and baseline covariates have been proposed. In an article
published in this journal, we develop an adaptation of parametric
G-computation.5 This involves fitting a parametric outcome model to patient-
level data for one of the studies to predict and contrast potential outcomes
in another study. When assumptions hold, parametric outcome modeling offers
greater statistical precision and efficiency than MAIC.5 Nevertheless, this
does not imply that the former is inherently superior to the latter.

Where there is limited overlap between covariate distributions across studies,
outcome modeling methods extrapolate the association between covariates and
outcome. In an insightful commentary, Vo warns of the dangers of extrapolating
a misspecified outcome model.6 Valid model-based extrapolation depends on
strong assumptions and relies on accurately capturing the true covariate-
outcome relationships. Vo advocates for the use of MAIC over G-computation
because the former does not extrapolate, and outlines the pitfalls of using
data-adaptive methods to estimate the outcome model.6

We congratulate Vo on his excellent editorial and thank the Editor for
arranging this exchange. In this counter-response, we examine four important
considerations for the development of covariate adjustment methodologies in
the context of indirect treatment comparisons. Firstly, we consider potential
advantages of weighting versus outcome modeling, placing focus on bias-
robustness. Secondly, we outline why model-based extrapolation may be required
and useful, in the specific context of indirect treatment comparisons with
limited overlap. Thirdly, we describe challenges for covariate adjustment
based on data-adaptive outcome modeling. Finally, we offer further
perspectives on the promise of doubly-robust covariate adjustment frameworks.

##  1 Bias-robustness: the case for weighting

Insofar, covariate-adjusted indirect comparison methods have relied on
parametric assumptions for unbiased estimation. In practice, these assumptions
will fail to some extent. Therefore, we must think about the bias-robustness
of different approaches. In parametric G-computation, bias arising from model
misspecification cannot be easily verified in the data. As highlighted by Vo
in his first simulated example, an outcome model that seems approximately
correct when fitted to the patient-level data would not necessarily fit well
in extrapolated regions.6 In practice, candidate models are compared using
diagnostic plots, information criteria or predictive performance measures such
as the mean squared prediction error. These measures cannot detect model
misspecification in extrapolated regions. Most model performance measures
average over the observed covariate distribution, and are not necessarily
transportable between different covariate distributions.7

In weighting approaches, evaluating the reduction in bias due to covariate
adjustment is more straightforward. Weighting methods do not necessarily
require knowledge of the trial assignment process. They seek weights that
attain cross-study balance in the distributions of measured covariates. By
checking for adequate balance after weighting, one can assess whether bias due
to differences in observed covariates has been mitigated.

Vo suggests adopting the individual-level covariate generation step proposed
by Remiro-Azócar et al.5 to fit the trial assignment model, typically a
logistic regression, using maximum-likelihood estimation.6 This approach
explicitly models the conditional probability of trial assignment as a
function of the baseline covariates, with odds weights derived from the
estimated probabilities.8 The typical model-building strategy is outcome-blind
and evaluates different specifications on a trial-and-error basis. One may
begin with a model containing additive terms for the relevant covariates, and
consider including non-linear transformations or higher-order terms, e.g.
polynomials, interactions, if imbalances remain.9

MAIC follows a different paradigm for weighting. It directly enforces balance
in covariate moments without explicitly modeling the conditional probability
of trial assignment. Covariate balance is viewed as a convex optimization
problem and the weights are solutions to such problem. MAIC is an entropy
balancing technique.10, 11 Entropy balancing approaches are more stable,
efficient and robust to model specification than the standard “inverse
weighting” modeling approaches.10, 12

Recent results imply that MAIC, balancing for covariate means, is
asymptotically unbiased if covariates have a linear relationship with the mean
difference (in the anchored scenario) or the absolute outcome means (in the
unanchored scenario), even if the implicit trial assignment model is
incorrectly specified.10, 12 Nevertheless, MAIC may be subject to bias if
there are non-linearities that are not accounted for by the, also implicit,
outcome model. MAIC remains biased in Vo’s first simulated example, partly
because there is a squared covariate-by-treatment product term in the outcome-
generating model, and only first-order moments (sample means) and not second-
order moments (sample variances) are balanced for.6 Because second-order
balance is enforced by balancing the means of the squared covariates,1
balancing both first- and second-order moments would have provided some
protection against bias. However, there is a risk that MAIC is unable to find
a solution to the convex optimization problem, and there may also be a
substantial cost in terms of precision.13

In summary, model misspecification bias is not easily diagnosed for outcome
modelling, particularly when there is extrapolation. Conversely, balance
diagnostics for weighting offer some indication on the level of bias reduction
due to covariate adjustment. It is important to make a distinction between
modeling and balancing approaches to weighting. Generally, the latter are less
susceptible to bias and offer greater stability and precision than the former.

##  2 Model-based extrapolation: undesirable but necessary

Where overlap is poor, weighting is inefficient. Effective sample sizes are
low and results may be sensitive to the inordinate influence of a few extreme
weights. A drop in precision is expected for all covariate adjustment methods
as overlap decreases. In our original simulation study, such reduction is more
marked for MAIC than for parametric G-computation.5 When modeling assumptions
hold, the latter is more precise and efficient than the former across all
scenarios.5 The increase in precision is achieved by extrapolating into non-
overlapping regions of the covariate space. The extrapolation is made using
strong parametric assumptions, which may be unfeasible. In our original
simulation study, variance estimates and coverage rates are valid for
G-computation with a correctly specified outcome model.5 Even so, the interval
estimates do not typically reflect the extrapolation uncertainty. This could
be problematic when modeling assumptions fail, particularly given that the
appropriate characterization of uncertainty is a central objective of HTA.14

Vo argues that weighting approaches reflect uncertainty more honestly.6 This
point is also raised in the discussion of our article,5 following Vansteelandt
and Keiding.15 The degree of extrapolation made by parametric G-computation is
implicit; it cannot be controlled and may hide any underlying lack of overlap.
Conversely, the presence of extreme weights in MAIC explicitly manifests the
high estimation uncertainty in data-sparse regions.

Vo recommends reducing the impact of extreme weights without resorting to
extrapolation.6 A suggested approach is weight truncation. Remiro-Azócar has
explored truncation in combination with MAIC, and shown that it can produce
substantial precision and efficiency improvements when sample sizes are low
and the extremity of the weights is high.16 However, truncation induces
asymptotic bias and is prone to ad hoc heuristics, e.g. deciding on cutoff
thresholds.16 Given the arbitrary nature of such decisions, weight truncation
is arguably not “honest” in quantifying uncertainty.

Weight stabilization is also recommended by Vo.6 The version of MAIC in our
original article fits a weighted marginal model of outcome on a time-fixed
binary treatment.5 Because this model is “saturated” — there is a parameter
for each unique predictor combination —- stabilized and unstabilized weights
yield identical results.17 Jackson et al. have developed an adaptation to MAIC
that minimizes the dispersion of the weights, improving precision but inducing
some bias.18 Remiro-Azócar recently proposes a modular extension to MAIC,
which improves precision and efficiency without resorting to bias-variance
trade-offs, particularly when sample sizes are low and covariates are highly
prognostic.16 It does so by estimating the treatment assignment mechanism in
the study with patient-level data, and is therefore not applicable in the
unanchored case.

The aforementioned approaches are useful for variance reduction, yet depend on
weighting solutions being available in the first place. Where there is no
common support, feasible weighting solutions providing satisfactory balance
may not exist.19 If the observed range of a covariate in the study with
patient-level data does not cover the respective mean published for the
comparator study, MAIC suffers from convergence failures and cannot even
produce an estimate.18 According to Vo, “we might need to accept that the
difference between populations is too large to be adjustable”.6 Nevertheless,
treatment effect estimates are needed for HTA and, in many jurisdictions, to
populate the economic models required to make reimbursement decisions. Vo
views model-based extrapolation as a flaw.6 In our view, extrapolation beyond
the overlapping regions of the covariate space may be necessary when the study
samples are divergent.

Finally, in the anchored scenario, Vo suggests excluding covariates that
overlap poorly but are not effect modifiers from the MAIC adjustment,6 as
indicated by current guidance.2 Our simulation study assumes a best-case
situation where these covariates can be identified. However, this is
challenging in practice, due to the scale dependence of effect modification
and to limited knowledge on the drivers of treatment effect heterogeneity,
particularly for novel therapies. As highlighted by Vo in his second simulated
example, excluding weak effect modifiers in order to reduce standard errors
may induce bias.6 There is a tension between satisfying the conditional
constancy or transportability of effects across studies and maintaining
overlap. Accounting for a greater number of covariates increases the
plausibility of the former but decreases the likelihood of the latter.

In conclusion, parametric G-computation relies on model-based extrapolation to
overcome limited overlap and improve precision with respect to MAIC. Model-
based extrapolation depends on strong assumptions and is typically viewed as
an undesirable feature. Nevertheless, it may be necessary in practice to
generate the treatment effect estimates that are required for HTA decision-
making.

##  3 The perils of data-adaptive G-computation

For asymptotically unbiased estimation, parametric G-computation relies on the
correct specification of the outcome model. Parametric models impose
restrictive assumptions on functional forms; for instance, that effects are
linear and additive on some transformation of the conditional outcome mean.
These assumptions may not be plausible where there are multiple continuous
covariates and complex non-linear relationships, particularly if background
theory is weak. Making strong parametric assumptions has a cost. If the
parametric model is incorrectly specified, the G-computation estimator is
subject to bias and this bias does not decrease with sample size, at any rate.
This may lead to substantial undercoverage: as the sample size grows, the
probability that the interval estimates contain the target estimand shrinks to
zero.

To reduce the risk of bias resulting from parametric model misspecification,
one may consider using data-adaptive methods, e.g. non-parametric techniques
or machine learning, to estimate the conditional outcome expectations.20 These
approaches make weaker modeling assumptions, but the flexibility also comes at
a cost. From a frequentist perspective, that is a slower rate of convergence
than the parametric rate.21 A misspecified parametric estimator attains fast
$n^{1/2}$ convergence, but to the wrong target value. Data-adaptive estimators
are more likely to converge to the true expectations, but converge more
slowly. This may lead to larger-than-desirable bias and poor coverage in
finite samples.

Consider the naive use of lasso-based G-computation in the setting explored by
Vo.6 The lasso is a regularization estimator that trades off variance for bias
to avoid overfitting and improve out-of-sample predictive performance. This
results in slow convergence rates, with bias and mean square error diminishing
slowly with sample size, particularly in high-dimensional settings. In Vo’s
second simulated example,6 the lasso dismisses various covariates that are
weak effect modifiers but are correlated with study assignment, ultimately
leading to bias.

In addition to slower than parametric convergence rates, data-adaptive
G-computation estimators also have non-normal sampling distributions.22 When
using parametric outcome models, the central limit theorem and the law of
large numbers can be applied to obtain convenient asymptotic approximations.
This enables one to readily construct variance and interval estimates, e.g.
through resampling methods such as the bootstrap and Wald normal interval
equations. Conversely, the asymptotic distribution of data-adaptive estimators
often is irregular and cannot be derived. This makes uncertainty
quantification difficult. Even where one can establish the asymptotic
distribution, the use of the bootstrap to compute interval estimates is not
supported by theoretical results as it requires certain convergence rate
conditions to hold.23

We propose moving away from a frequentist framework and exploring the use of
Bayesian G-computation. Non-parametric outcome models such as Bayesian
additive regression trees (BART)24, 25, 26 are very flexible and well-suited
to capture complex functional forms. Because BART directly samples from the
posterior distribution of the treatment effect, principled variance and
interval estimation are possible without having to resort to the bootstrap.
Bayesian methods naturally integrate the analysis into a probabilistic
formulation that is desirable for HTA,14 and allow the incorporation of
substantive prior information.

BART is implemented in many software packages, requires little parameter
tuning, and has performed excellently in causal inference applications.25, 27
Nevertheless, BART is not a panacea. As per the lasso, it regularizes heavily
to penalize overfitting, which may induce non-negligible bias into the
estimation of treatment effects.26 Little is known theoretically about BART’s
asymptotic (frequentist) properties. In the causal inference literature, the
method has exhibited poor coverage in scenarios with limited overlap because
it does not extrapolate well.28 Regression trees are piecewise step functions
that partition the covariate space into disjoint cuts. Consequently, they
extrapolate beyond the overlap region using a constant value.28

In conclusion, the blind use of data-adaptive estimation together with
G-computation may lead to bias in finite samples and is constrained by limited
theoretical justification for valid statistical inference. Data-adaptive
G-computation estimators are more flexible than their parametric counterparts.
Nevertheless, flexible estimators still extrapolate where overlap is limited.
The extrapolations of very flexible models may be intrinsically flawed,
impairing point and interval estimation under poor overlap.

##  4 The appeal of doubly-robust adjustment

Vo proposes the use of doubly-robust approaches such as augmented weighting
estimators,6 in conjunction with the individual-level covariate simulation
step by Remiro-Azócar et al.5 Augmented weighting methods apply two working
models: a model for trial assignment conditional on covariates, and another
for the outcome conditional on treatment and covariates. Only one of the two
models needs to be correctly specified to achieve asymptotically unbiased
estimation. In general, doubly-robust estimators should be less prone to model
misspecification bias than singly-robust estimators; they offer two
opportunities for valid adjustment instead of one, mitigating the risks of
applying the wrong outcome model. Nevertheless, the two working models are
typically parametric. In practice, there is seldom adequate background
knowledge on the true trial assignment and outcome-generating processes.
Therefore, both parametric models may be incorrect, in which case doubly-
robust adjustment is subject to bias.

A valuable suggestion by Vo is the application of data-adaptive methods within
a doubly-robust framework.6 An oft-ignored advantage of such framework is that
it allows for the slower convergence of working models.29, 30 Augmented
weighting approaches combine two singly-robust estimators, such that the
overall convergence rate is as fast or faster than the convergence rate of
each separate estimator. Consider that a flexible non-parametric regression or
machine learning technique is used to fit both working models. Even if each
working model converges at a rate slower than the optimal parametric $n^{1/2}$
rate, the product of the individual rates may still be $n^{1/2}$ or faster.29,
30 Cross-validation based ensembles of data-adaptive techniques (e.g. Super
Learner) also satisfy the rate condition if one of the candidate estimators
does.31

Valid inference for data-adaptive doubly-robust estimators remains a
challenge. Their asymptotic normality, which allows for simple approaches to
construct variance and interval estimates, relies on an unverifiable “Donsker”
condition from empirical process theory.32 Such condition may be inappropriate
in high-dimensional settings, restricting the complexity of the data-adaptive
methods that are permitted.

The Donsker condition can be relaxed using sample-splitting procedures such as
cross-fitting.29, 33 These create an additional source of variation in finite
samples, and need to be iterated repeatedly to remove seed-dependence.29, 33
This is computationally burdensome if there is no closed-form expression for
the variance, as bootstrapping would incur considerable run-time. The
asymptotic distribution of estimators can be learned by characterizing their
efficient influence function, as suggested by Vo.6 Nevertheless, we suspect
that finding tractable approximations is not trivial with limited access to
patient-level data.

In summary, doubly-robust approaches require the correct specification of
either a trial assignment or an outcome model to obtain consistent estimation.
While singly-robust data-adaptive estimators may suffer from slow convergence
rates, doubly-robust estimators can accommodate the use of slower-converging
data-adaptive models. Nevertheless, their asymptotic normality relies on an
unverifiable Donsker condition, which restricts the complexity of the models
that are permitted. Such condition could be overcome through the use of
sample-splitting, but this may require substantial computational resources
when closed-form expressions for the variance are unavailable.

##  5 Conclusions

Weighting is the most popular approach for covariate adjustment in the context
of indirect treatment comparisons with limited patient-level data. Recently,
alternatives based on outcome modeling have been developed. Such methods
enable the estimation of treatment effects in “limited overlap” settings where
MAIC cannot be applied. Under no failures of assumptions, parametric outcome
modeling increases precision and efficiency with respect to weighting,
particularly when overlap is poor. It does so by extrapolating outside the
overlap region, which involves untestable and potentially untenable parametric
assumptions. While functional forms may hold approximately in the overlap
region, this is not necessarily the case when extrapolating beyond it.

Specifying an appropriate outcome model – one that seems approximately correct
in the overlap region but that also extrapolates well – is challenging in many
settings, and can quickly become an overwhelming process. Analysts are
typically drawn to simple but restrictive models, assuming linear effects to
be the norm. Subject-matter knowledge about the outcome-generating process is
often fragile. Data-adaptive outcome modeling techniques are attractive to
mitigate the risk of bias arising from parametric model misspecification.
These may not require the pre-specification of functional forms and can
provide some automation to the model selection process. Nevertheless, the slow
convergence rates of flexible models may also lead to bias. In addition,
flexible estimators can extrapolate poorly and may offer limited theoretical
justification for valid uncertainty propagation.

Precision and efficiency under no failures of assumptions should not be the
only characteristics determining the choice of a suitable estimator. In
general, the level of bias reduction due to covariate adjustment is easier to
assess for weighting than for outcome modeling. Several authors argue that
entropy balancing approaches to weighting such as MAIC are more bias-robust
than outcome modeling methods.6, 10, 12 Doubly-robust methods that combine a
trial assignment model with a model for the conditional outcome expectation
offer two opportunities for correct model specification, and are likely more
bias-robust than both MAIC and parametric G-computation. The doubly-robust
framework also allows for the use of slower-converging data-adaptive models,
although adequate variance and interval estimation may have its challenges.

Standard (unadjusted) indirect treatment comparisons rely on a very strong
assumption: the unconditional constancy or transportability of relative
effects across studies. Covariate-adjusted indirect comparisons seek to relax
this assumption by adjusting for baseline covariates. In doing so, they rely
on the correct specification of parametric models. The selection of an
appropriate model is challenging and subject to many “researcher degrees of
freedom”. The development of innovative covariate adjustment methods that help
to reduce these degrees of freedom is imperative. The valid characterization
of uncertainty, a key requirement of HTA,14 is another important factor in
determining the suitability of these methods for decision-making. The
performance of novel covariate adjustment methods should be examined in
simulation studies reflecting typical scenarios in the context of indirect
treatment comparisons with limited patient-level data.

Finally, innovative methodologies may not overcome aspects such as limited
overlap and the unavailability (or inconsistent definition) of influential
covariates. The former may require model-based extrapolation, which always
calls for a leap of faith. The latter can be addressed through the development
of core patient characteristic sets that define important prognostic
covariates to be measured and reported among specific therapeutic areas.34

## Acknowledgments

The authors thank Tat-Thang Vo for his insightful editorial and the editors of
Research Synthesis Methods for facilitating this discussion and inviting us to
reply. Anna Heath is funded by Canada Research Chair in Statistical Trial
Design; Natural Sciences and Engineering Research Council of Canada (award No.
RGPIN-2021-03366).

### Financial disclosure

Funding agreements ensure the authors’ independence in writing and publishing
the article.

### Conflict of interest

Antonio Remiro-Azócar is employed by Bayer plc but declares no conflicts of
interest. The views expressed in this article do not necessarily represent
those of Bayer plc. Anna Heath and Gianluca Baio declare no potential
conflicts of interest.

### Data Availability Statement

Data sharing is not applicable to this article as no datasets were generated
or analyzed.

## References

  * 1 Phillippo D, Ades T, Dias S, Palmer S, Abrams KR, Welton N. NICE DSU technical support document 18: methods for population-adjusted indirect comparisons in submissions to NICE. 2016.
  * 2 Phillippo DM, Ades AE, Dias S, Palmer S, Abrams KR, Welton NJ. Methods for population-adjusted indirect comparisons in health technology appraisal. Medical Decision Making 2018; 38(2): 200–211.
  * 3 Signorovitch JE, Wu EQ, Andrew PY, et al. Comparative effectiveness without head-to-head trials. Pharmacoeconomics 2010; 28(10): 935–945.
  * 4 Phillippo DM, Dias S, Elsada A, Ades A, Welton NJ. Population adjustment methods for indirect comparisons: a review of national institute for health and care excellence technology appraisals. International journal of technology assessment in health care 2019; 35(3): 221–228.
  * 5 Remiro-Azócar A, Heath A, Baio G. Parametric G-computation for Compatible Indirect Treatment Comparisons with Limited Individual Patient Data. Research synthesis methods 2022; 13(6): 716–744.
  * 6 Vo TT. A cautionary note on the use of G-computation in population adjustment. Research synthesis methods 2023; 14(3): 338–341.
  * 7 Steingrimsson JA, Gatsonis C, Dahabreh IJ. Transporting a prediction model for use in a new target population. In press, American Journal of Epidemiology 2023; 192(2): 296–304.
  * 8 Westreich D, Edwards JK, Lesko CR, Stuart E, Cole SR. Transportability of trial results using inverse odds of sampling weights. American journal of epidemiology 2017; 186(8): 1010–1014.
  * 9 Kyle RP, Moodie EE, Klein MB, Abrahamowicz M. Evaluating flexible modeling of continuous covariates in inverse-weighted estimators. American journal of epidemiology 2019; 188(6): 1181–1191.
  * 10 Josey KP, Berkowitz SA, Ghosh D, Raghavan S. Transporting experimental results with entropy balancing. Statistics in medicine 2021; 40(19): 4310–4326.
  * 11 Phillippo DM, Dias S, Ades A, Welton NJ. Equivalence of entropy balancing and the method of moments for matching-adjusted indirect comparison. Research synthesis methods 2020; 11(4): 568–572.
  * 12 Cheng D, Tchetgen ET, Signorovitch J. On the Double-Robustness and Semiparametric Efficiency of Matching-Adjusted Indirect Comparisons. In press, Research Synthesis Methods 2022\.
  * 13 Campbell H, Park JE, Jansen JP, Cope S. Standardization allows for efficient unbiased estimation in observational studies and in indirect treatment comparisons: A comprehensive simulation study. arXiv preprint arXiv:2301.09661 2023\.
  * 14 Claxton K, Sculpher M, McCabe C, et al. Probabilistic sensitivity analysis for NICE technology assessment: not an optional extra. Health economics 2005; 14(4): 339–347.
  * 15 Vansteelandt S, Keiding N. Invited commentary: G-computation–lost in translation?. American journal of epidemiology 2011; 173(7): 739–742.
  * 16 Remiro-Azócar A. Two-stage matching-adjusted indirect comparison. BMC Medical Research Methodology 2022; 22(1): 1–16.
  * 17 Shiba K, Kawahara T. Using propensity scores for causal inference: pitfalls and tips. Journal of epidemiology 2021; 31(8): 457–-463.
  * 18 Jackson D, Rhodes K, Ouwens M. Alternative weighting schemes when performing matching-adjusted indirect comparisons. Research Synthesis Methods 2021; 12(3): 333–346.
  * 19 Phillippo DM, Dias S, Ades A, Welton NJ. Assessing the performance of population adjustment methods for anchored indirect comparisons: A simulation study. Statistics in Medicine 2020; 39(30): 4885–4911.
  * 20 Le Borgne F, Chatton A, Léger M, Lenain R, Foucher Y. G-computation and machine learning for estimating the causal effects of binary exposure statuses on binary outcomes. Scientific reports 2021; 11(1): 1–12.
  * 21 Naimi AI, Mishler AE, Kennedy EH. Challenges in obtaining valid causal effect estimates with machine learning algorithms. In press, American Journal of Epidemiology 2021\.
  * 22 Vaart v. dA. Higher order tangent spaces and influence functions. Statistical Science 2014; 29(4): 679–686.
  * 23 Bickel PJ, Götze F, Zwet vWR. Resampling fewer than n observations: gains, losses, and remedies for losses. Statistica Sinica 1997; 7(1): 1–31.
  * 24 Chipman HA, George EI, McCulloch RE. Bayesian ensemble learning. Advances in neural information processing systems 2007; 19: 265.
  * 25 Hill JL. Bayesian nonparametric modeling for causal inference. Journal of Computational and Graphical Statistics 2011; 20(1): 217–240.
  * 26 Hahn PR, Murray JS, Carvalho CM. Bayesian regression tree models for causal inference: Regularization, confounding, and heterogeneous effects (with discussion). Bayesian Analysis 2020; 15(3): 965–1056.
  * 27 Kern HL, Stuart EA, Hill J, Green DP. Assessing methods for generalizing experimental impact estimates to target populations. Journal of research on educational effectiveness 2016; 9(1): 103–127.
  * 28 Zhu AY, Mitra N, Roy J. Addressing positivity violations in causal effect estimation using Gaussian process priors. Statistics in Medicine 2023; 42(1): 33–51.
  * 29 Chernozhukov V, Chetverikov D, Demirer M, et al. Double/debiased machine learning for treatment and structural parameters. The Econometrics Journal 2018; 21(1): C1–-C68.
  * 30 Zeng Z, Kennedy EH, Bodnar LM, Naimi AI. Efficient Generalization and Transportation. arXiv preprint arXiv:2302.00092 2023\.
  * 31 Vaart AWvd, Dudoit S, Laan MJvd. Oracle inequalities for multi-fold cross validation. Statistics & Decisions 2006; 24(3): 351–371.
  * 32 Van Der Vaart AW, Wellner JA, Vaart v. dAW, Wellner JA. Weak convergence.  Springer .  1996\.
  * 33 Zivich PN, Breskin A. Machine learning for causal inference: on the use of cross-fit estimators. Epidemiology (Cambridge, Mass.) 2021; 32(3): 393–401.
  * 34 Vo TT, Vuong ML, Tu PHT, Duong KL. Development of minimum reporting sets of patient characteristics in epidemiological research: a methodological systematic review. medRxiv 2023: 2023–02.


# On cubic Thue equations and the common index divisors of cyclic cubic fields

Mohammed Seddik
Université d’Évry Val d’Essonne, Université Paris-Saclay
Laboratoire de Math matiques et Modélisation d’Évry (UMR 8071)
I.B.G.B.I., 23 Bd. de France, 91037 Évry Cedex, France
<EMAIL_ADDRESS>
###### Abstract

In this paper, we investigate the common index divisors of cyclic cubic
fields. Let $a,b,c,d$ and $k$ are integers, we then solve the following Thue
cubic equations:

$ax^{3}+bx^{2}y+cxy^{2}+dy^{3}=k\ $

when $a,bc+d$ are odd and $3$ doesn’t divide $v_{2}(k)$.

††2010 _Mathematics Subject Classification_ : 11R04; 11DXX; 11R16; 11R33;
11R09.†† _Key words and phrases_ : Cubic Thue equations, cyclic cubic fields,
common index divisors of cyclic cubic fields.

##  1 Introduction and statement of main result

In 1909, Thue [26] proved the following important result : If
$f(x,y)\in\mathbb{Z}[x,y]$ is an irreducible binary form of degree at least 3,
and $k$ non-zero integer, then the Diophantine equation $f(x,y)=k$ has only a
finite number of integer solutions $(x,y)$. However, the method used does not
help in finding the values of $x$ and $y$ in question. It is only possible to
obtain an upper bound for the number of solutions. Baker [5] has shown that
bounds can also be given for the magnitude of the solutions. The bounds for
$\arrowvert x\arrowvert$ and $\arrowvert y\arrowvert$ are extremely big, even
for small values of degree of $f(x,y)$.

If $f(x,y)$ is an irreducible binary cubic form with negative discriminant,
Delauney [9] and Nagell [23] showed that the equation $f(x,y)=1$ has at most
five integer solutions $(x,y)$. Now if its discriminant is positive, then
Evertse [11] showed that the equation $f(x,y)=1$ has at most twelve integer
solutions $(x,y)$. Recently, Bennett [6] refined Delauney-Nagell-Evertse
result as follows: if $f(x,1)$ has at least two distinct complex roots, then
the equation $f(x,y)=1$ possesses at most $10$ solutions in integers $x$ and
$y$.

Again, let $f(x,y)$ be an irreducible binary cubic form, the general equation
$f(x,y)=k$ with $|D(f)|>\gamma k^{33}$, $k$ non-negative integer, and $\gamma$
certain positive constant, Siegel proved there is at most $18$ solutions if
$D(f)$ is positive (resp. Fjellstedt proved there is at most $14$ solutions
when $D(f)$ is negative). See [23, p.208].

In 1984, Ayad [3] proved that if $f(x,y)$ is a binary form of degree $3$ with
coefficients in $\mathbb{Z}$, $Aut(f)$ its automorphisms group and $H(f)$ its
Hessian, then $Aut(f)$ is trivial except when $H(f)=\lambda g(x,y)$,
$\lambda\in\mathbb{Z}^{*}$ and $g(x,y)$ is equivalent to $x^{2}-xy-y^{2}$. In
this last case, $Aut(f)$ is cyclic of order $3$ and $f$ is equivalent to one
binary form of type :

$f_{m,n}(x,y)=mx^{3}-nx^{2}y-(n+3m)xy^{2}-my^{3},\quad m,n\in\mathbb{Z},$

so, the number of representations of integer $k$ by $f(x,y)$ is divisible by
$3$. Note that the case $m=k=1$ is proved by Avanesov [2].

Recently, in [29] Wakabayashi proved that for $k=1$, the Thue equation
$f_{m,n}(x,y)=1$ has at most three integer solutions except for a few known
cases. Using the Padé approximation method he obtained a lower bound for the
size of solutions.

However, only few results are known when the group of automorphisms is
trivial. The family of cubic Thue equations $f_{1,n}(x,y)=\pm 1$, with
$n\geq-1$, was studied by Thomas in [25]. Using Baker’s method, Thomas proved
that it has only trivial solutions except for a finite number of values of the
parameter $n$, explicitly for $n<10^{8}$. Later, Mignotte solved the remaing
cases for this family of equations in [19].

Let $n$ be a rational integer and $\mathbb{K}_{n}=\mathbb{Q}(\theta)$ be a
cyclic cubic number field generated by a root $\theta$ of
$f_{1,n}(X,1)=X^{3}-nX^{2}-(n+3)X-1$ and let $\mathbb{O}_{\mathbb{K}_{n}}$ be
its ring of integers. The polynomial $f_{1,n}(X,1)$ has discriminant
$(n^{2}+3n+9)^{2}$. If $n^{2}+3n+9$ is square-free, then we have the
discriminant of $\mathbb{K}_{n}$, $D(\mathbb{K}_{n})=(n^{2}+3n+9)^{2}$ and
$\mathbb{O}_{\mathbb{K}_{n}}=\mathbb{Z}[\theta]$ (there exist infinitely many
such $n$, cf. Cusick [8, pp. 63-73]).

For $\theta$ a root of $P_{n}$, Pëtho and Lemmermeyer in [16] proved that for
all $\alpha\in\mathbb{Z}[\theta]$ either $|N(\alpha)|\geq 2n+3$, or $\alpha$
is associated to an integer. Moreover, if $|N(\alpha)|=2n+3$, then $\alpha$ is
associated to one of the conjugates of $\alpha-1$. In [20], Mignotte, Pëtho
and Lemmermeyer, by Baker’s method and used the results of [16] and [8,
pp.63-73] solved completely the case $m=1$, $n\geq-1$ and $1\leq k\leq 2n+3$.
Lettl, Pëtho and Voutier [17], following an idea of Chudnovsky [7], improved
the usual estimate of the Padé approximation and gave a good upper bound for
the size of solutions for $f_{1,n}(x,y)=k$, $n\geq 30$ and $k$ is arbitrary.
Wakabayashi [30] studied Thue inequality $|f_{m,n}(x,y)|\leq k$ with two
parameters $m,n$.

In 2011, A. Hoshi [13] studied the case when $k$ is a positive divisor of
$n^{2}+3n+9$, and gave a correspondence between integer solutions to the
parametric family of cubic Thue equations

$x^{3}-nx^{2}y-(n+3)xy^{2}-y^{3}=k$

and isomorphism classes of the simplest cubic fields. For more details on the
study of simplest cubic fields see [24].

A. Togbé [28], using Baker’s method and the results obtained by L. C.
Washington [31] and O. Lecacheux [15], solved the family of parametrized Thue
equations

$x^{3}-(n^{3}-2n^{2}+3n-3)x^{2}y-n^{2}xy^{2}-y^{3}=\pm
1,\quad\text{when}\;n\geq 1.$

A. Togbé [27] using Baker’s method and the results obtained by Y. Kishi [14],
solved the family of parametrized Thue equations

$x^{3}-n(n^{2}+n+3)(n^{2}+2)x^{2}y-(n^{3}+2n^{2}+3n+3)xy^{2}-y^{3}=\pm
1,\quad\text{when}\;n\geq 0.$

Recently, Balady [4] showed that there are many families of cubics polynomials
irreducibles with square discriminants. We consider the family of polynomials

$X^{3}+(n^{7}+2n^{6}+3n^{5}-n^{4}-3n^{3}-3n^{2}+3n+3)X^{2}+(-n^{4}+3n)X-1$

used in [4] which satisfy the hypothesis of our main result.

Our method of proof is different to the previous methods, we use the common
index divisors of cyclic cubic fields. Now we state our main result.

######  Theorem 1.

Let $a,b,c,d$ and $k$ be integers where $3\nmid v_{2}(k)$ and $a,bc+d$ are
odd. Suppose that the polynomial $aX^{3}+bX^{2}+cX+d$ is irreducible with
square discriminant. Then the cubic Thue equations

$ax^{3}+bx^{2}y+cxy^{2}+dy^{3}=k$

has no integer solution $(x,y)$.

Below we mention some immediate consequences of our result on the resolution
of a parametric families of cubic equations.

######  Corollary 1.

Let $k\in\mathbb{Z}$ where $3\nmid v_{2}(k)$. Then the parametric families of
cubic Thue equations

  1. 1.

$mx^{3}-nx^{2}y-(n+3m)xy^{2}-my^{3}=k,\;m,n\in\mathbb{Z},\;m\;\text{odd},\;\gcd(m,n)=1,$

  2. 2.

$x^{3}-(n^{3}-2n^{2}+3n-3)x^{2}y-n^{2}xy^{2}-y^{3}=k,\;\;n\in\mathbb{Z},\;n\not\equiv
1\;mod\;4,$

  3. 3.

$x^{3}-n(n^{2}+n+3)(n^{2}+2)x^{2}y-(n^{3}+2n^{2}+3n+3)xy^{2}-y^{3}=k,\;n\in\mathbb{Z},$

  4. 4.

$x^{3}+(n^{7}+2n^{6}+3n^{5}-n^{4}-3n^{3}-3n^{2}+3n+3)x^{2}y+(-n^{4}+3n)xy^{2}-y^{3}=k,\;n\in\mathbb{Z},$

have no integer solution $(x,y)$.

######  Corollary 2.

Let $a,b,c,d$ and $k$ be integers where $3\nmid v_{2}(k)$ and $a,bc+d$ are
odd. Suppose that the polynomial $aX^{3}+bX^{2}+cX+d$ is irreducible with
square discriminant. Then the homogeneous form

$ax^{3}+bx^{2}y+cxy^{2}+dy^{3}=kz^{3}$

has only integer solution $(x,y,z)=(0,0,0)$.

######  Corollary 3.

Let $a,b,c,d$ and $t$ be integers where $v_{2}(t)=1$ and $a,bc+d$ are odd.
Suppose that the polynomial $aX^{3}+bX^{2}+cX+d+t$ is irreducible with square
discriminant. Then the equation

(1) $\displaystyle ax^{3}+bx^{2}+cx+d=ty^{2}$

has no integer solution $(x,y)$.

The equation (1) was studied by Mordell [21, pp.255-261]. He proved the
following important result: if $a,b,c,d,t$ are arbitrary integers and the
polynomial $ax^{3}+bx^{2}+cx+d$ has no squared linear factor in $x$, then the
equation (1) has only a finite number of integer solutions. Here in under the
by hypothesis of Corollary 3, the equation has no integer solution.

## References

  * [1] E. T. Avanesov, _A bound for the number of representations by a special class of binary cubic forms of positive discriminant_ , Acta Arith. 20 (1972), 17-31.
  * [2] , _Solution in integers of the indeterminate equation $x^{3}-2x^{2}y-5xy^{2}-y^{3}=1$_, Ivanov. Gos. Ped. Inst. Ucen. Zap. 61 (1969), 61-83.
  * [3] M. Ayad, _Automorphismes d’une forme binaire cubique et représentation d’entiers_ , C. R. Acad. Sci. Paris Sér. I Math. 299 (1984), no. 19, 959-962.
  * [4] S. Balady, _Families of cyclic cubic fields,_ J. Number Theory. 167 (2016), 394-406.
  * [5] A. Baker, _Contribution to the theory of Diophantine equations,_ Phil. Trans. Roy. Soc. London, Ser. A, 167 (1967/68), 717-208.
  * [6] M. A. Benett, _On the representation of unity by binary cubic forms_ , Trans. American Math. soc. Vol. 353 (2000), N∘ 4, 1507-1534.
  * [7] G. V. Chudnovsky, _On the Thue-Siegel,_ Ann. of Math. 117 (1983), 325-382.
  * [8] T. W. Cusick, _Lower bound for regulators_ , in Number Theory Noordwijkerhout Springer Verlag (1983), 63-73, Ed : H. Jager, LNM Vol. 1068, Springer Verlag (1984).
  * [9] B. Delauney, _Über die Darstellung der Zahlen durch die binären kubischen Formen von negative Discriminante_ , Math. Z. 31 (1930), 1-26.
  * [10] V. Ennola, _Cubic number fields with exceptional units_ , Computational number theory (Debrecen, 1989), de Gruyter, Berlin, 103-128.
  * [11] J. H. Evertse, _On the Representation of Integers by Binary Cubic Forms of Positive Discriminant_ , Invent. Math. 73 (1983), 117-138.
  * [12] C. Heuberger, A. Pëtho and R. F. Tichy, _Complete solution of parametrized Thue equations_ , Acta Math. Inform. Univ. Ostraviensis. 6 (1998), 91-113.
  * [13] A. Hoshi, _On correspondence between solutions of a family of cubic Thue equations and isomorphism classes of the simplest cubic fields_ , J. Number Theory. 131 (2011), 2135-2150.
  * [14] Y. Kishi, _A family of cyclic polynomials whose roots are systems of fundamental units_ , J. Number Theory 102  (2003), 90-106.
  * [15] O. Lecacheux, _Units in number fields and elliptic curves_ , Advances in Number Theory, Oxford Science Publications, Oxford University Press, New York (1993), 293-301.
  * [16] F. Lemmermeyer, A. Pëtho, _Simplest cubic fields_ , Manuscripta Math. 88 (1995), 285-321.
  * [17] G. Lettl, A. Pëtho, P. Voutier, _Simple families of Thue inequalities_ , Trans. Amer. Math. Soc. 351 (1999), 1871-1894.
  * [18] P. Llorente, E. Nart, _Effective determination of the rational primes in a cubic field_ , Proc. Amer. Math. Soc. 87 (1983), 579-585.
  * [19] M. Mignotte, _Verification of a conjecture of E. Thomas_ , J. Number Theory 44 (1993), 172-177.
  * [20] M. Mignotte, A. Pëtho, F. Lemmermeyer, _On the family of Thue equations $x^{3}-(n-1)x^{2}y-(n+2)xy^{2}-y^{3}=k$_, Acta Arith. 76 (1996), no. 3, 245-269.
  * [21] L. J. Mordell, _Diophantine Equations_ , Academic Press, London, (1969).
  * [22] P. Morton, _Characterizing cyclic cubic extensions by automorphism polynomials,_ J. Number Theory. 49 (1994), 183-208.
  * [23] T. Nagell, _Remarques sur une classe d’équations indéterminées,_ Arkiv för Math. 8 (1969), 199-214.
  * [24] D. Shanks, _The simplest cubic fields_ , Math. Comp. 28 (1974), 1137-1152.
  * [25] E. Thomas, _Complete solutions to a family of cubic diophantine equations_ , J. Number Theory. 34 (1990), 235-250.
  * [26] A. Thue, _Über Annäherungswerte algebraischer Zahlen_ , J. reine angew. Math. 135 (1909), 284-305.
  * [27] A. Togbé, _On the solutions of a parametric family of cubic Thue equations_ , Bull. Braz. Math. Soc. 39 (4) (2008), 537-554.
  * [28] , _A parametric family of cubic Thue equations_ , J. Number Theory. 107 (2004), 63-79.
  * [29] I. Wakabayashi, _Number of solutions for cubic Thue equations with automorphisms_ , Ramanujan J. 14 (2007), 131-154.
  * [30] , _Simple families of Thue inequalities_ , Ann. Sci. Math. Comp. 31 (2007), 211-232.
  * [31] I. Washington, _A family of cubic fields and zeros of 3-adic L-functions_ , J. Number Theory. 63 (1997), 408-417.