TITLE: How to re-estimate a theoretical outcome knowing rng is imperfect
QUESTION [0 upvotes]: I have 3 regular dice with six sides. I am playing a game where I win when I match all three dice to the same side. My expected shannon bits of entropy H(X) evaluates to 7.754887502 bits. However, my dice aren't very good as a random number generator and the entropy is actually measured at 6.965784285. How do I recalculate my new theoretical win rate?
$$
P_{expected} = (\frac 1 6) ^ 2
$$

REPLY [0 votes]: First, I believe $P_{expected}=(\frac{1}{6})^2$ because it doesn't matter what the first die is, just that the other 2 match. The value you give for $P_{expected}$ is the value for a specific value, say, 6. There are six such values, so $(\frac{1}{6})^3•6=(\frac{1}{6})^2$.Second, there is a difference between theoretical probability and experimental probability. I wouldn't recalculate my theoretical probability; I would re-run my experiment again (and again and again), updating the experimental probability by averaging the results. If your experimental probability continues to differ significantly from your theoretical probability, I wouldn't change the theoretical model; I would get better dice. :)


TITLE: Finding the kernel of $\alpha: K[X,Y,Z]^{3}\rightarrow \langle X,Y,Z\rangle$, $(f,g,h)\mapsto Xf+Yg+Zh$.
QUESTION [4 upvotes]: I am trying to do exercise $2.3$ of Reid's "Undergraduate Commutative Algebra":

Let $A=K[X,Y,Z]$ where $K$ is a field, and $m=\langle X,Y,Z\rangle$. I have to show that the kernel of the homomorphism
  $$
\alpha: A^{3}\rightarrow m, (f,g,h)\mapsto Xf+Yg+Zh
$$
   is generated by three elements of $A^{3}$.

I think the elements are $(-y,x,0), (-z,0,x)$ and $(0,-z,y)$. It is clear that
$\langle(-y,x,0), (-z,0,x),(0,-z,y)\rangle\subseteq\ker(\alpha)$, but I am not able to prove the other inclusion.

REPLY [4 votes]: We know $Xf+Yg+Zh=0$. By using the long division algorithm one can write $g=Xg_1+g_2$, and $h=Xh_1+h_2$ with $g_2,h_2\in K[Y,Z]$. Plugging this into $Xf+Yg+Zh=0$ we get $Xf+XYg_1+Yg_2+XZh_1+Zh_2=0$. This shows that $Yg_2+Zh_2\in (X)$, and since $g_2,h_2$ don't contain $X$ we get $Yg_2+Zh_2=0$. Thus $Y\mid h_2$ and $Z\mid g_2$, and then $g_2=-Zu$, $h_2=Yu$, with $u\in K[Y,Z]$. Then $f+Yg_1+Zh_1=0$. Now one can write $(f,g,h)=(-Y,X,0)g_1+(-Z,0,X)h_1+(0,-Z,Y)u$.


TITLE: Disjunctive normal form (BOTH dnf and cnf) example help
QUESTION [1 upvotes]: Consider
$$T \vee ( \neg x \wedge y ) \vee ( x \wedge y )\tag{1}$$
In class we went over examples of how many things can be both DNF and CNF... e.g.,

$(\neg z) \vee y$

can be thought of as both

$(\neg z) \vee y\qquad\quad\,$ DNF
$((\neg z) \vee y) \wedge T\quad$ CNF

Can someone explain why example (1) cannot work like that, i.e., as in being one large conjunction of a disjunction as in (2)?
$$T \wedge (T \vee ( \neg x \wedge y ) \vee ( x \wedge y ))\tag{2}$$
Is it due to having more than one disjunction inside?

REPLY [1 votes]: For DNF :

a formula is a disjunctive normal form if and only if it is a disjunction of one or more conjunctions of one or more literals.

Tus, the formula :

$(\lnot z) \lor y$

is in DNF beacuse is the disjunction of two "degenerate" conjunctions : $\lnot z$ and $y$.
For CNF, we have that :

a formula is a conjunctive normal form (CNF) if and only if it is a conjunction of clauses, where a clause is a disjunction of literals.

Thus :

$((\lnot z) \lor y) \land T$

is in CNF, because we have the disjunction of the two literals : $\lnot z$ and $y$ that makes the clause : $(\lnot z) \lor y$; then a second clause : $T$ (a "degenerate" disjunction) and the two disjunctions are conjoined into : $((\lnot z) \lor y) \land T$.
Regarding the formula :

$T \land [T \lor (\lnot x \land y ) \lor ( x \land y )]$

we have three conjunctions : $(\lnot x \land y )$ and $( x \land y )$ and the "degenerate" : $T$, and they are disjoined into :

$T \lor (\lnot x \land y ) \lor ( x \land y )$;

this formula is in DNF; but adding the part : $T \land \ldots$, the resulting formula is nor more a DNF.

With the formula :

$T \land [T \lor (\lnot x \land y ) \lor ( x \land y )]$

we can apply Distributivity to its subformula:
$$[(\lnot x \land y ) \lor ( x \land y )] \equiv [y \land (x \lor \lnot x)] \equiv (y \land T) \equiv y$$
to get the equivalent :

$[T \land (T \lor y)] \equiv (T \land T) \equiv T$

which is both a CNF and a DNF equivalent to the original formula.


text:
\begin{document}

\begin{abstract}
  Modern combustion engines incorporate a number of actuators and sensors
  that can be used to control and optimize the performance and
  emissions. We describe a semi-automatic method to simultaneously measure
  and calibrate the actuator settings and the resulting behavior of the
  engine. The method includes an adaptive process for refining the
  measurements, a data cleaning step, and an optimization procedure. The
  optimization works in a discretized space and incorporates the conditions
  to describe the dependence between the actuators and the engine behavior
  as well as emission bounds. We demonstrate our method on
  practical examples.
\end{abstract}

\maketitle

\section{Introduction}
\label{sec:intro}
Due to the wish to save fossil fuels, stringent maximal emission limits and challenging customer preferences, modern internal combustion engines (ICEs) become more and more complex.
Indeed, research towards the improved construction of combustion engines is highly relevant for reaching the \CT emission targets set by the European Union, cf.\ \cite{TERM2016}.

This engine development results in an increasing number of \emph{actuators} and \emph{sensors}.
There are currently in the order of ten different actuators and sensors each.
Examples for actuators include the amount of injected fuel, exhaust recirculation control, air valve angle, etc.
Sensors measure, e.g., the temperature, maximal point of the cylinder, torque, exhaust emission, etc.
The actuators allow to produce a certain torque and revolution frequency, which describe the two main requirements on the engine in usage.
However, several different settings of the actuators can result in the same
torque/revolution frequency combination.
Moreover, their dependence is involved and not known exactly a priori, i.e., it has to be measured and approximated.
In this paper we deal with the \emph{optimal engine calibration problem}, i.e., to efficiently approximate this dependence by few measurements and to choose optimal actuators settings.

The engine calibration problem consists of determining a so-called
\emph{engine manifold}, which determines for each torque/revolution frequency
combination a corresponding setting for actuators. This manifold is usually
discretized and the resulting \emph{solution map} is hard-coded into the engine control unit. The settings on the engine manifold
have to be chosen in such a way that they are consistent across different
torque/frequency combinations, i.e., the engine manifold should be continuous.
Moreover, the resulting settings need to obey several
restrictions in order to avoid damage of the engine as well as bounds on the emissions produced.
The emissions are measured with respect to so-called \emph{driving cycles}.
These are certain prescribed changes in the torque/frequency settings over time which should resemble usage in practice.
These driving cycles are applied to the engine on a measurement bench and
the resulting emissions have to be bounded.

The calibration process described above involves two main steps.
First, the dependence between the actuator settings and the output has to be measured.
A naive enumeration of a uniform grid of possible actuator settings and
interpolation would require exponential time in the number of actuators.
Therefore, one needs to design a process to perform the measurements at interesting areas in order to speed up the measurement process.
Additionally, the actuator settings are continuously changed without waiting until a steady-state has been reached.
This allows to save time, but also introduces measurement
errors like hysteresis effects.
Second, one should optimally choose the resulting engine manifold and solution map.
Thus, based on the measurements, optimal actuator
settings obeying the above-mentioned restrictions should be chosen in order
to obtain an approximation of an engine manifold.

In this article we propose a way to solve the engine calibration problem,
which consists of the following contributions:
\begin{myitemize}
\item \emph{Adaptive meshing:} The density of measurements is adapted
  within areas where the measured function is sensitive with respect to its
  inputs, while keeping the density of measurements coarse where it does
  not. This leads to a more accurate representation of the engine manifold than
  with a uniform grid approach with a fraction of the measurements.

\item \emph{Data cleaning:} Before optimization, the measured data are
  cleaned by filtering out redundancies and noise.

\item \emph{Discretization and Optimization:}
We discretize the space of measurands in a fashion that fits the format of lookup tables as they are stored in the engine control unit. Using discrete optimization techniques, we select among the measured (and cleaned) data such actuator settings that minimize fuel consumption of the engine while its pollutant emissions conform to current regulations. The selected settings are drivable in the sense that the actuator's variation speed is bounded in order to prevent engine damage.
\end{myitemize}

We would like to stress that the algorithm described in this article is a \emph{tool for}
the engineer, not a replacement. While it runs automatically once its
parameters have been set, the setting of these parameters, e.g., the
determination of the subset of actuators which is to be varied in a given
situation, requires extensive knowledge and experience. In this sense it is
a semi-automatic process.
\smallskip

This article is structured as follows.
In Section~\ref{sec:state-of-the-art} we provide a brief review of the relevant literature.
A high-level description of the mathematical problem is given in Section~\ref{sec:math}.
From this, we arrive at the corresponding steps of our process.
Section~\ref{sec:calibration} provides the details of our method.
In Section~\ref{sec:avl} we present a practical case-study, and Section~\ref{sec:experiments} contains the experimental results.
We close with a few concluding remarks.

\section{State of the Art}\label{sec:state-of-the-art}

In this section we describe the state of the art of calibration methods for the optimization of ICEs. Many commercial and research products exist for the measurement and calibration of ICEs.
In the following we give a brief overview of the different measurement and optimization approaches and the products that implement them.

\subsection{Determining Data Boundaries}

Throughout its operation the measurand values of an ICE have to stay within certain boundaries.
Some of these boundaries are set in place to avoid destructive situations, e.g., for cylinder pressures or critical device temperatures.
Others are induced, e.g., by emission and noise regulations.
The set of actuator settings such that neither of the boundary constraints is violated is called the operational space. It can be detected by different techniques.

One classical method is the measurement via star like measurement paths with a defined step width. Each measurement path or ray will be interrupted if one observable violates a destructive or nondestructive limit. The convex hull of all measurement path points then defines the operational space. This technique is typically automatized and part of most commercial tools, like \Cameo \cite{Gschweitl2001} and \TopExpert \cite{TOPexpert2004}.

The \MATLAB toolbox \LOLIMOT (local linear model tree) contains a further method, developed by TU Darmstadt and Porsche AG. This tool determines the boundary limits iteratively. The measurement starts with a coarse detection of the operational space. In the following iteration steps, the algorithm defines linear measurement paths, which are orthogonal to the facets of the convex hull of the measured boundary points. These measurement paths will be interrupted if hard or soft boundary conditions are violated \cite{Schueler2000}.

Besides the determination of the operational space, other programs directly
monitor the boundary conditions. For example, BMW developed the tool \mbminimize, which uses regression models for engine limits \cite{Knoedler2003}. IAV created \RapidMeasurement, which facilitates faster measurements by employing mathematical and physical models to compensate latency effects \cite{Roepke2007}.

\subsection{Test Run Planning }
\label{Full Factorial}

Generating a \textit{uniform grid} of measurements is the most obvious method to study the behavior of an unknown system.
However, the cost of this approach is exponential in the number of actuators.

The most common technique for test run planning is the \emph{Design of Experiments (DoE)} approach.
The DoE is a test run planning strategy which employs mathematical modeling to describe the behavior of a system of observables with a minimal number of measurement points.
The classical test plans are grouped into central composite, space-filling, and certain designs which are called \emph{D-optimal}.
The latter approach is followed, e.g., by \LOLIMOT, which generates local models via linear neuronal networks \cite{Isermann2010}.
The tool \mbminimize combines a D-optimal measurement plan with additional interactive features \cite{Knoedler2003}.

\subsection{Measurement Methods}

The measurement methods can be divided into stationary, quasi-stationary and transient/dynamic measurement. The stationary measurement method is most frequently used. Dynamic measurement techniques are still limited to very small application areas, e.g., the dynamic optimization of emission objectives \cite{Schueler2000}.

During a quasi-stationary measurement, the parameters are varied slowly. Therefore the whole physical system follows the dynamic parameters with a bounded contouring error. The data acquisition permanently records the behavior of the whole system. By varying one or more parameters, one observes the impact instantaneously. By inverting the measurement path it is possible to compensate the contouring error \cite{Isermann2010}. A local model derivation is possible for up to three parameters \cite{Isermann2010}.

\DYNMET is a measurement strategy, in order to train dynamic neuronal networks. The system is stimulated with a variation on every input signal. For nonlinear processes one uses APRBS (amplitude modulated pseudo random binary signal) \cite{Isermann2010}.

\subsection{Model-based Optimization}

One approach to the calibration process is to approximate the engine behavior by a mathematical model and then to optimize this model exactly.
However, ICEs are difficult to describe by mathematical or physical models that cover the whole operational space.
Especially models for emission predictions are strongly limited.
ICEs are systems which follow the \emph{Multiple Input--Multiple Output (MIMO)} paradigm.
Models often do not contain physical correlations, but describe the engine as a black box system \cite{Isermann2010}.

Typically, one uses polynomial functions and \emph{least squares} techniques in order to approximate a simple model to measured data.
Yet strongly nonlinear regions are often not described well enough by these approaches \cite{Mitterer}.

A further approach are artificial neuronal networks (ANN), which work well for nonlinear regimes. Furthermore, this black-box approach does not need special knowledge about the ICE \cite{Mitterer}. Neuronal networks are separated mainly into
Multi-Layer-Perceptrons (MLPs) and radial base function (RBF) networks. The principle of \LOLIMOT is similar to RBF networks \cite{Isermann2010}.

The toolbox \ASCMO uses the superposition of weighted base functions \cite{KRUSE2010739}.
The weights are derived automatically by a statistical learning approach (SLA).
This method is robust and needs only a small number of data points to train the network.
The toolbox also includes a calibration suite for the so-derived models.

As a key drawback model-based approaches, be it via neural networks or polynomial approximations, often fail to recognize critical behavior.
This is due to the fact that the models usually rely on smooth curves which are fitted to a given set of data, whereas the processes underlying the combustion are usually nonsmooth.
Hence, for engines calibrated via model-based optimization, the actual fuel consumption and emission values may substantially deviate from the predictions.
This is one key motivation for the present work.

\section{Mathematical Problem Description}
\label{sec:math}

\noindent
The calibration of an internal combustion engine is the procedure to derive
an optimal engine manifold. To reduce the complexity of the model and to
enhance practical implementation, these manifolds are discretized to obtain
solution maps (see Section~\ref{sec:discretization} below).

In our setting, knowledge of the engine behavior with respect to variations of its actuator settings is obtained by means of physical experiments on a test bench. In addition to the revolution frequency, typical actuators include the injected fuel quantity, the injection angle, or the valve pressure.
The generated torque of the engine is a measurand.

Technically there is a relevant distinction between direct and controlled actuators. Direct actuators
such as injected fuel quantity, injection angle, or valve pressure can be set directly on the engine, while controlled actuators are set indirectly. For example, the revolution frequency is a controlled actuator that is regulated via a brake on the engine shaft. However, in our mathematical model this distinction is not relevant. Here all actuator setting are considered as coordinates of points in the admissible domain $\domain$, while measured values such as the engine torque and exhaust emissions yield coordinates of points in the admissible range $\codomain$.

The aforementioned side constraints include limits on pollutant emission as well as physical requirements such as engine temperature limits.
They necessitate that not only torque, but several other output values of the engine are measured as well.
A realistic engine model features $m\geq 8$ actuators and $n\geq 14$ measurands.
In our mathematical model we represent the relation between the setting of
$m$ actuators and $n$ sensor values by a function
\[
F\colon \R^m \to \R^n.
\]
Throughout we make the fundamental assumption that $F$ is continuous, but not necessarily differentiable everywhere.
Further, we will assume the actuator values to be restricted to a box $\domain\subset \R^m$, which we call the \emph{admissible domain}.
The noncritical sensor values define another box $\codomain\subset \R^n$, the \emph{admissible range}.
The exact definitions of $\domain$ and $\codomain$ will be stated in the subsequent Section \ref{sec:init}.
The \emph{feasible space} of $F$ is the set
\[
\feasible \define \{ (u,y)\in \domain\times\codomain \suchthat y=F(u) \} \ \subseteq \ \R^m\times\R^n.
\]
The actuators correspond to the coordinates of $u$; examples are the revolution frequency and the amount of fuel injected.
Typical sensor values, i.e., coordinates of $F(u)$, include the torque and the emission of carbon monoxide.

We denote by ``$\frequency$'' the index of the actuator for revolution frequency and by ``$\torque$'' the index of the torque sensor values.
Then we call
\[
\OP \define \{(u_\frequency,y_\torque) \suchthat (u,y)\in\feasible\}
\]
the \emph{operation field} of $F$.
The operation field represents the behavior of the engine with respect to revolution frequency and torque.
While our methods are more general, we focus on this particular pair of actuator and sensor values in our analysis.

\begin{figure}
	\centering

	\definecolor{gray1}{gray}{0.87}
	\definecolor{gray2}{gray}{0.6}
	\definecolor{gray3}{gray}{0.75}
	\begin{tikzpicture}[scale=1]

	\draw [<->,thick] (0,6.1) node (yaxis) [above] {$\torque$}
	|- (7.5,0) node (xaxis) [right] {$\frequency$};

	\draw (1,0) -- (1,5.6);
	\draw (2,0) -- (2,5.6);
	\draw (3,0) -- (3,5.6);
	\draw (4,0) -- (4,5.6);
	\draw (5,0) -- (5,5.6);
	\draw (6,0) -- (6,5.6);
	\draw (7,0) -- (7,5.6);

	\draw (0.5,0) -- (0.5,5.6);
	\draw (1.5,0) -- (1.5,5.6);
	\draw (2.5,0) -- (2.5,5.6);
	\draw (3.5,0) -- (3.5,5.6);
	\draw (4.5,0) -- (4.5,5.6);
	\draw (5.5,0) -- (5.5,5.6);
	\draw (6.5,0) -- (6.5,5.6);

	\draw (0,0.8) -- (7,0.8);
	\draw (0,1.6) -- (7,1.6);
	\draw (0,2.4) -- (7,2.4);
	\draw (0,3.2) -- (7,3.2);
	\draw (0,4) -- (7,4);
	\draw (0,4.8) -- (7,4.8);
	\draw (0,5.6) -- (7,5.6);
	\draw (0,0.4) -- (7,0.4);
	\draw (0,1.2) -- (7,1.2);
	\draw (0,2) -- (7,2);
	\draw (0,2.8) -- (7,2.8);
	\draw (0,3.6) -- (7,3.6);
	\draw (0,4.4) -- (7,4.4);
	\draw (0,5.2) -- (7,5.2);

	\draw [fill=gray1] (0.5,0.8) rectangle (1,1.2);
	\draw [fill=gray1] (1,0.8) rectangle (1.5,1.2);
	\draw [fill=gray1] (1.5,0.8) rectangle (2,1.2);
	\draw [fill=gray3] (2,0.8) rectangle (2.5,1.2);
	\draw [fill=gray2] (2.5,0.8) rectangle (3,1.2);
	\draw [fill=gray2] (3,0.8) rectangle (3.5,1.2);
	\draw [fill=gray2] (3.5,0.8) rectangle (4,1.2);
	\draw [fill=gray3] (4,0.8) rectangle (4.5,1.2);
	\draw [fill=gray2] (4.5,0.8) rectangle (5,1.2);
	\draw [fill=gray3] (5,0.8) rectangle (5.5,1.2);
	\draw [fill=gray1] (5.5,0.8) rectangle (6,1.2);
	\draw [fill=gray1] (6,0.8) rectangle (6.5,1.2);

	\draw [fill=gray1] (1,1.2) rectangle (1.5,1.6);
	\draw [fill=gray1] (3,1.2) rectangle (3.5,1.6);
	\draw [fill=gray2] (4,1.2) rectangle (4.5,1.6);

	\draw [fill=gray1] (1,1.6) rectangle (1.5,2);
	\draw [fill=gray1] (3,1.6) rectangle (3.5,2);
	\draw [fill=gray1] (4,1.6) rectangle (4.5,2);

	\draw [fill=gray1] (1,2) rectangle (1.5,2.4);
	\draw [fill=gray1] (1.5,2) rectangle (2,2.4);
	\draw [fill=gray1] (2,2) rectangle (2.5,2.4);
	\draw [fill=gray1] (2.5,2) rectangle (3,2.4);
	\draw [fill=gray1] (3,2) rectangle (3.5,2.4);
	\draw [fill=gray1] (3.5,2) rectangle (4,2.4);
	\draw [fill=gray1] (4,2) rectangle (4.5,2.4);

	\draw [fill=gray1] (1,2.4) rectangle (1.5,2.8);
	\draw [fill=gray1] (2,2.4) rectangle (2.5,2.8);
	\draw [fill=gray1] (3,2.4) rectangle (3.5,2.8);
	\draw [fill=gray1] (4.5,2.4) rectangle (5,2.8);

	\draw [fill=gray1] (1,2.8) rectangle (1.5,3.2);
	\draw [fill=gray1] (2,2.8) rectangle (2.5,3.2);
	\draw [fill=gray1] (3,2.8) rectangle (3.5,3.2);
	\draw [fill=gray1] (4.5,2.8) rectangle (5,3.2);

	\draw [fill=gray1] (1.5,3.2) rectangle (2,3.6);
	\draw [fill=gray1] (2,3.2) rectangle (2.5,3.6);
	\draw [fill=gray3] (2.5,3.2) rectangle (3,3.6);
	\draw [fill=gray1] (3,3.2) rectangle (3.5,3.6);
	\draw [fill=gray3] (3.5,3.2) rectangle (4,3.6);
	\draw [fill=gray3] (4,3.2) rectangle (4.5,3.6);
	\draw [fill=gray1] (4.5,3.2) rectangle (5,3.6);
	\draw [fill=gray1] (5,3.2) rectangle (5.5,3.6);

	\draw [fill=gray1] (2,3.6) rectangle (2.5,4);
	\draw [fill=gray1] (5.5,3.6) rectangle (6,4);

	\draw [fill=gray1] (2.5,4) rectangle (3,4.4);
	\draw [fill=gray1] (3,4) rectangle (3.5,4.4);
	\draw [fill=gray1] (3.5,4) rectangle (4,4.4);
	\draw [fill=gray1] (4,4) rectangle (4.5,4.4);
	\draw [fill=gray1] (4.5,4) rectangle (5,4.4);
	\draw [fill=gray1] (5,4) rectangle (5.5,4.4);
	\draw [fill=gray1] (5.5,4) rectangle (6,4.4);
	\draw [fill=gray1] (6,4) rectangle (6.5,4.4);

	\draw [fill=gray1] (4,4.4) rectangle (4.5,4.8);
	\draw [fill=gray1] (4.5,4.4) rectangle (5,4.8);
	\draw [fill=gray1] (5,4.4) rectangle (5.5,4.8);

	\end{tikzpicture}
	\caption{\label{fig:frog} Schematic of discretized operation field and driving cycle. Resistance times are represented by gray shades. (Also compare Figures \ref{fig:nox-nedc} and \ref{fig:nox-random}.)}
\end{figure}

\subsection{Continuous Optimization Problem}\label{sec:continuous-opt}

For the sake of a concise exposition we now describe an idealized
continuous optimization problem which has the actual optimization
problem that we want to solve as a natural discretization.

\subsubsection*{Optimization Space}

The feasible region of our continuous optimization problem is given by all engine manifolds.
Each such manifold is given as the image of a map $M \colon \OP \to \R^m$ that assigns actuator settings to a given frequency and torque value pair in the operation field.
A basic requirement is that these maps are continuous.

Moreover there are vital additional conditions to consider.
Varying actuators too fast might damage the engine.
Therefore in the (continuous) final solution map the variation speed of every actuator is bounded by constants $\Delta_\actuator$ for $\actuator \in [m]$, where $[m] \define \{1, \dots,
m\}$.
This yields the \textit{drivability constraint} for a given map $M$ with respect to the actuator~$\actuator$:
\begin{align}\label{eq:driv-constraint}
  \abs{M(f, t)_{\actuator} - M(f', t')_{\actuator}} \leq \Delta_\actuator \cdot \norm{(f, t),(f', t')},
\end{align}
for all $(f, t)$, $(f', t') \in \OP$, where $\norm{\cdot}$ is some norm.
Thus, the map $M_\actuator:\OP\to\R$ is Lipschitz continuous with constant $\Delta_\actuator$.
We define $\Omega$ as the set of all maps~$M$ which are feasible in the sense that they obey the drivability constraint with respect to each actuator.

\subsubsection*{Driving Cycle and Emission Constraints}

In accordance with current government regulations and common test cycles, the engine behavior is optimized with respect to pre-defined scenarios, known as driving cycles.
In our continuous model a \emph{driving cycle} is a time-parametrized curve
\[
\dcCont \colon [0,1] \to \OP,
\]
whose purpose is to simulate the phases of acceleration and constant speed of real-world driving patterns.
For the sake of the simplicity of the exposition here we will assume that a fixed driving cycle is given; cf.\ Figure~\ref{fig:rnd-profile} for an example.
Yet it is also natural to take the combination of several driving cycles into account, and our approach also covers this slightly more general situation.

An important constraint prescribed by regulations is to bound the resulting emissions along the driving cycle.
Emission pollutants include carbon monoxide (\CO), hydrocarbons (\HC), nitrogen oxides (\NX) as well as
particulate matter (\PM) and number (\PN). Denote by $\emission$ the index
set that corresponds to pollutant emissions and by $\emconstraint_p$ the
emission limit for pollutant~$p$ over the driving cycle.

Consider an engine manifold map $M \in \Omega$ and a pollutant~$p \in \emission$.
The frequency and torque values along $\dcCont$ produce actuator settings $M(\dcCont(t))$, which result in output values $F(M(\dcCont(t)))$, including the pollutant~$p$.
The integral of these values must satisfy the \textit{emission constraint} for all $p \in \emission$:
\begin{align}\label{eq:em-constraint}
  \int_{0}^{1} [F(M(\dcCont(t)))]_p \,\text{d}t \leq \emconstraint_p.
\end{align}
For instance, Table \ref{tab:euro-norms} shows the diesel engine emission constraints for \euro norms 3-6c (E3-E6c).

\subsubsection*{The optimization formulation}

Using the terminology developed so far, the continuous optimization problem is
\begin{align}\label{eq:continuous-opt}
  \min_{M \in \Omega}\ &  \int_{0}^{1} [\dcCont(t)]_\fuel \,\text{d}t\\
  & \int_{0}^{1} [F(M(\dcCont(t)))]_p \,\text{d}t \leq \emconstraint_p &
  \forall p \in \emission.\notag
\end{align}
The optimization goal is to minimize the total consumed fuel along the given driving cycle~$\dcCont$.

\renewcommand{\arraystretch}{1.3}
\begin{table}
	\caption{\label{tab:euro-norms} Selection of \euro emission constraints for passenger cars with compression ignition engine.}
	\small
	\centering
	\begin{tabular*}{\textwidth}{@{}l@{\extracolsep{\fill}}lrrrrr@{}}\toprule
		Em. & Unit &  E3 (2001) &  E4 (2006) &  E5a (2011) & E6b (2015) & E6c (2018) \\ \midrule
		\CO & $\text{mg}\,/\,\text{km}$ & 2300 & 1000 & 1000 & 1000 & 1000 \\
		\HC & $\text{mg}\,/\,\text{km}$ & 200 & 100 & 100 & 100 & 100 \\
		\NX & $\text{mg}\,/\,\text{km}$ & 150 & 80 & 60 & 60 & 60 \\
		\PN & $\text{mg}\,/\,\text{km}$ & -- & -- & 5 & 4.5 & 4.5 \\
		\PM & $$1\,/\,\text{km}$$ & -- & -- & -- & $6\cdot 10^{12}$ & $6\cdot 10^{11}$ 	\\\bottomrule
	\end{tabular*}
\end{table}

\subsection{Discretization} \label{sec:discretization}

To obtain a finite-dimensional problem, the model is discretized to yield
\emph{(characteristic) engine maps}. We consider combinations of~$k$
revolution frequencies and $k$ torque demands and subdivide the operation
field $\OP$ of~$F$ into $k^2$ congruent rectangles $\OP_{ft}$, where $f$,
$t \in[k]$ denote the rectangle's frequency and torque coordinate,
respectively. Thus, engine maps are given as $k {\times} k$-matrices which
for each combination yield the corresponding value of a particular
actuator.  They are stored permanently in the engine control unit.
A \emph{solution map} consists of a complete set of engine maps, one for each actuator.
In this way a solution map yields a discretization of the continuous solution map described in Section~\ref{sec:continuous-opt}.

Then, in our terminology, the solution map takes as input a frequency and
torque pair $(f,t)\in [k]\times[k]$ and yields as output an admissible
actuator setting.
The latter is a point $u\in\domain$ with $(u,F(u))\in\feasible$
such that the pair $(u_\frequency, y_\torque)$ lies in the rectangle of the
discretized operation field corresponding to the input coordinates $(f,t)$.

Our goal is to obtain a solution map which is optimal with respect to a
given objective, while conforming to several constraints.  Since typical
constraints are continuous but nonlinear and the solution map itself is
discrete, this optimization problem belongs to the wide class of
mixed-integer/discrete nonlinear optimization problems, which are often
difficult to handle.

Our goal can be formulated as follows: Find a finite set
\[
\dmeasurements \define \{(u^1,y^1),(u^2,y^2),\dots,(u^N,y^N)\}
\]
of $N$ points in $\feasible$, i.e., a set of points $u^q\in\domain$ such that $y^q=F(u^q)\in\codomain$, which allows to compute a sufficiently good approximation to an optimal solution $M$ of the continuous optimization problem \eqref{eq:continuous-opt}.
The points $(u^q,y^q)$ are obtained via actual measurements, and thus $\dmeasurements$ is called the \emph{data set}.
Elements of $\dmeasurements$ are called \textit{data points}.
The final selection of the data points is the result of a cycle of measurements and optimization steps, which are the discussed in the sections below.

We will assume that none of the points in the data set lies on the boundary of any rectangle $\OP_{ft}$.
In this case, the discretization of the operation field into $k^2$ rectangles partitions $\dmeasurements$ into sets
\[
  S_{ft} = \bigl\{ (u,y) \in \dmeasurements \suchthat (u_{\frequency},y_{\torque}) \in\OP_{ft} \bigr\}.
\]
We call each set $S_{ft}$ a \emph{stack}, and the entire partition $\kOP \define \{S_{ft} \suchthat f,t\in[k]\}$ is the \emph{$k$-operation field} of $F$ with respect to $\dmeasurements$.
In the solution map each entire rectangle $\OP_{ft}$ will be represented by a single measurement $d_{ft}\in S_{ft}$.
We will postpone the question how to pick these representatives.

\subsubsection*{Discrete drivability constraint}

Due to the uniform discretization of the $k$-operation field, \textit{a discrete drivability constraint} merely has to bound the difference between actuator settings of representatives of neighboring rectangles $\OP_{ft}$. That is
\begin{align}\label{eq:driv-constraint-disc}
	\abs{[d_{ft}]_\actuator - [d_{gs}]_\actuator} \leq \Delta_\actuator \qquad \text{for all } \actuator\in[m],
\end{align}
and all tuples $(f,t)$, $(g,s)$, where either $g = f \pm 1$ or $s = t\pm 1$.

\subsubsection*{Discrete emission constraint}

To each rectangle $\OP_{ft}$ we further associate a weight $\omega_{ft}$.
This weight is set to zero if the intersection of $\OP_{ft}$ with the image of the curve $\dcCont$ is empty.
Otherwise $\omega_{ft}$ is a positive value that reflects the resistance time, i.e., the duration of the curve $\dcCont$ staying in the rectangle $\OP_{ft}$.
The function
\[
  \dcDisc \colon [k]\times[k]\to\R_{\geq 0},\ (f,t)\mapsto \omega_{ft}
\]
serves as a discrete analogue of the continuous driving cycle $\dcCont$.
Note that the map $\dcDisc$ only records which parts of the operating field are met by  $\dcCont$, but it ignores the order in which this happens.

The \textit{discrete emission constraint} is now given as
\begin{equation}\label{eq:em-constraint-disc}
\sum_{(f,t)\in[k]\times[k]} \omega_{ft}\cdot [d_{ft}]_p \ \leq \emconstraint_p \qquad\text{for all } p\in\emission.
\end{equation}

The practical driving cycles, for instance the New European Driving Cycle (\nedc) or Real World Driving Cycle (\rand), are given by operational profiles which map a time to a gear-velocity combination; cf.\ Figure \ref{fig:rnd-profile}.
The weights $\omega_{ft}$ can be derived from these profiles. The revolution frequency at a certain time can be calculated directly from the current gear-velocity combination. For the requested engine torque not only the speed but also the acceleration has to be taken into account, along with the car's mass, roll drag and air flow resistance.
\begin{figure}[t]
	\includegraphics[width=0.48\textwidth]{Random_velocity.pdf}
	\hfill
	\includegraphics[width=0.48\textwidth]{Random_gear.pdf}
	\caption{Operational profile of \rand driving cycle.}
	\label{fig:rnd-profile}
\end{figure}

\subsubsection*{Discrete Optimization Problem}

In the discrete version of problem \eqref{eq:continuous-opt} we are given
the set $\dmeasurements$, subdivided into $k^2$ stacks $S_{ft}$ with
corresponding weights~$\omega_{ft}$. The optimization objective is now to pick a set of $k^2$ representatives $d_{ft}$ of the stacks $S_{ft}$ that solve
\begin{equation}\label{eq:discrete-opt}
  \min_d\ \sum_{(f,t)\in[k]\times[k]} \omega_{ft}\cdot [d_{ft}]_\fuel,
\end{equation}
subject to the discrete drivability constraint \eqref{eq:driv-constraint-disc} and the discrete emission constraint \eqref{eq:em-constraint-disc}.

\section{Semi-automatic Calibration}
\label{sec:calibration}
\noindent
The purpose of this section is to explain the algorithmic core of our method.
As its input it is given the domain $\domain$, equipped with a grid, and an evaluation oracle for the function $F$.
It returns a set $\SOL\subset\feasible$ of representatives of the $k$-operation field (recall that its cardinality equals $k^2$).
Algorithm~\ref{algo:calibration} is a first rough sketch, which will be detailed in the following.
In practice the evaluation oracle for $F$ is given by an engine mounted on a test-bench.

\begin{algorithm}
  \SetKwInOut{Input}{Input}
  \SetKwInOut{Output}{Output}

  \Input{\mbox{admissible domain $\domain$}, \mbox{grid $G$,} \mbox{admissible range $\codomain$,} \mbox{data set $\dmeasurements$,} \mbox{precision parameter $k$,} \mbox{evaluation oracle $F$}}
  \Output{solution map $\SOL$, updated data set $\dmeasurements$}

  initialization

  \While{solution map $\SOL$ incomplete}
  {
    iteration step: adds to $\dmeasurements$ and updates $G$

    data cleaning: reduces $\dmeasurements$

    select $\SOL$ from $\dmeasurements$ via integer linear program
  }
  \Return $(\SOL, \dmeasurements)$

  \caption{Engine calibration procedure}
  \label{algo:calibration}
\end{algorithm}

\subsection{Initialization}\label{sec:init}
The calibration procedure is initialized with the domain $\domain$ equipped with a grid $G$, the set of noncritical target values $\codomain$, a (possibly empty) set of data points $\dmeasurements$ and a precision parameter $k$.
The function $F$ is given implicitly; for a given setting of the actuators, the sensors yield the respective function value by means of a physical measurement.
We assume that the data points in $\dmeasurements$ reflect true values of $F$.

The domain
\begin{equation}\label{eq:domain-product}
  \domain \ =\ \prod_{i=1}^{m}I_i
\end{equation}
is a product of intervals $I_1,I_2,\dots,I_m$, where $I_i$ defines the range of variation of actuator $i$.
It may happen that $I_i$ degenerates to a single point.
Then the corresponding actuator is called \textit{static}, otherwise it is called \textit{dynamic}.

There are two cases to distinguish.
In the first case, the set $\dmeasurements$ of data points is empty.
Then, for each dynamic actuator $i$, the corresponding interval $I_i$ is subdivided into parts of equal lengths.
Otherwise, if $\dmeasurements$ is not empty, then we assume that each interval $I_i$ is equipped with its own subdivision.
In either case, the product form \eqref{eq:domain-product} induces a grid structure $G$ on the domain.
If $\dmeasurements$ is empty then this grid $G$ is uniform.
In the later stages of the optimization, however, $G$ will become more and more non-uniform.

Each measurand has an interval $J_j$ of noncritical values.
For example, the temperature of the test engine has to stay within certain bounds to prevent it from damage.
Then \[\codomain \ =\ \prod_{j=1}^{n}J_j\,.\]
In practice, $\codomain$ is given in part by the physical test engine (e.g., the aforementioned temperature limits) as well as by external factors such as government regulations (e.g., emission limits).
Just like $F$, which is given by the physical test engine, we will assume the set $\codomain$ to stay fixed throughout the whole calibration process.
The precision parameter $k$ is dictated by the engine control unit's engine map format.

During the various steps of the calibration, data points will be added to the set $\dmeasurements$.
Thus, it may happen that $\dmeasurements$ becomes prohibitively large to perform subsequent steps of the calibration.
How to weed out less relevant measurements is the subject of Section~\ref{subsec:cleaning} below.

\subsection{Iteration Step}\label{sec:iteration-step}

The basic iteration step can be subdivided into two phases: The generation of a measurement plan, followed by the actual measurement, which is combined with a refinement of the grid.

\subsubsection*{Generation of the Measurement Plan}

For our given grid $\grid$ and data set $\dmeasurements$, we construct an abstract graph $\graph=\graph(\grid, E)$ as follows.
The nodes of $\graph$ are the grid boxes determined by $\grid$.
Two $m$-dimensional grid boxes are joined by an edge if their intersection is a grid box of dimension $m-1$.
If there are no measurements yet, i.e., if $\dmeasurements$ is empty, then the grid $\grid$ is uniform, and the graph $\graph$ is the dual graph of a cubical cell complex.
Due to non-uniform refinement, the structure of $\graph$ will become more complicated.

The graph $\graph$ is equipped with nonnegative node and edge weights. For a grid box $B$, we denote by $\#B$ the number of points $(u^q,y^q)\in\dmeasurements$ such that $u^q\in B$. Then the weight of a grid box $B$ is chosen as
\begin{equation}\label{eq:wB}
  w(B) = \frac{\vol(B)}{\#B+1}.
\end{equation}
The quotient $\#B/\vol(B)$ is the \emph{data density} of $B$.
Similarly, we define the weight of an edge between two adjacent grid boxes $B$ and $B'$ as
\begin{equation}\label{eq:wBB}
  w(B,B') = \frac{\vol(B)+\vol(B')}{\#(B\cup B')+1}.
\end{equation}
There is some room for adjusting these weights; the general idea is that the weights should be roughly reciprocal to the data density.
But they should stay finite even if the data density is zero.

Due to the involved nature of the measurement process, it is not advisable to take arbitrary independent measurements.
In fact, the following concept is crucial as it reflects how the actual measurement process is organized on an engine test bench in practice.
\begin{definition}(Measurement Ramp)
  For two points $u^q$, $u^r\in\domain$, we call the set
  \[
  \Big\{u^q + i\cdot\frac{u^r-u^q}{\ell-1} \suchthat i=0,1,\dots, \ell-1\Big\}
  \]
  the \emph{measurement ramp} from $u^q$ to $u^r$ with $\ell$ measurements.
\end{definition}

Zero-entries of $u^r-u^q$ correspond to actuators which are locally static.
Suppose that there is an admissible point $u=u^M\in\domain$ where the last measurement took place.
If this does not exist, we choose $u$ uniformly at random in the domain $\domain$.
Let $B$ be the grid box containing~$u$.
We may assume that $B$ is unique, since $u$ has been constructed in a randomized fashion.
The two steps of the generation of the measurement plan are as follows:
\begin{enumerate}[I.]
\item\label{it:random-grid-box} \emph{Random grid box:}
  Pick a grid box $B'$ at random with probability $\tfrac{w(B')}{\overline w}$, where $\overline w$ is the sum of the weights of all boxes.
\item \emph{Measurement path:}
  Determine a shortest path $B_0=B,B_1,\dots,B_s=B'$ in $\graph$ from $B$ to $B'$ using Dijkstra's algorithm with respect to the weights; cf.\ \cite[\S2.2]{Cook1998}.
  In each box $B_q$ for $q \in [s]$, pick a point $u^{M+q}$ uniformly at random.
  Connect the points $u^M,u^{M+1},\dots,u^{M+s}$ by $s$ measurement ramps
  with~$\ell$ measurements each.
\end{enumerate}
This will result in a total of up to $s\cdot \ell$ new measurements to be added to the set $\dmeasurements$.

\begin{figure}
  \definecolor{gray1}{gray}{0.9}
  \definecolor{gray2}{gray}{0.7}
  \definecolor{gray3}{gray}{0.5}
  \definecolor{grayCross}{gray}{0.6}

  \begin{tikzpicture}[scale=1]

  \draw [<->,thick] (0,4.5) node (yaxis) [above] {Act. 2}
  |- (4.5,0) node (xaxis) [right] {Act. 1};

  \draw (1,0) -- (1,4);
  \draw (2,0) -- (2,4);
  \draw (3,0) -- (3,4);
  \draw (4,0) -- (4,4);
  \draw (0,1) -- (4,1);
  \draw (0,2) -- (4,2);
  \draw (0,3) -- (4,3);
  \draw (0,4) -- (4,4);
  \draw (1.5,1) -- (1.5,2);
  \draw (2.5,1) -- (2.5,3);
  \draw (3.5,2) -- (3.5,3);
  \draw (1,2.5) -- (2,2.5);
  \draw (2.5,2.5) -- (4,2.5);

  \draw [fill=gray1] (0,1) rectangle (1,2);
  \draw [fill=gray1] (3,0) rectangle (4,1);
  \draw [fill=gray1] (1,2) rectangle (2,2.5);
  \draw [fill=gray1] (2.5,2) rectangle (3,2.5);
  \draw [fill=gray2] (1,1) rectangle (1.5,2);
  \draw [fill=gray2] (2,1) rectangle (2.5,2);
  \draw [fill=gray3] (3,2.5) rectangle (3.5,3);
  \draw [fill=gray3] (1,2.5) rectangle (2,3);

  \draw (0.6,3.4) node[cross, line width=3pt,color=grayCross] {};
  \draw (3.7,1.8) node[cross, line width=3pt,color=grayCross] {};

  \draw[->, line width=1.2pt] (0.6,3.4) -- (0.4,2.3);
  \draw[->, line width=1.2pt] (0.4,2.3) -- (1.8,2.2);
  \draw[->, line width=1.2pt] (1.8,2.2) -- (1.9,1.2);
  \draw[->, line width=1.2pt] (1.9,1.2) -- (1.4,0.5);
  \draw[->, line width=1.2pt] (1.4,0.5) -- (2.8,0.4);
  \draw[->, line width=1.2pt] (2.8,0.4) -- (2.6,1.3);
  \draw[->, line width=1.2pt] (2.6,1.3) -- (3.7,1.8);

   \node[rectangle split, rectangle split parts=2, rectangle split part fill={gray3,gray1}, draw] at (4.5,4) {
   	\textcolor{white}{High Density}
   	\nodepart{two}
   	Low Density
   };
  \end{tikzpicture}
  \tikzstyle{every node}=[circle, draw, fill=black!50, inner sep=0pt, minimum width=4pt]
  \begin{tikzpicture}[scale=1.5]

  	\draw (0,0) -- (3,0);
  	\draw (0,1) -- (3,1);
  	\draw (0,3) -- (3,3);
  	\draw (0,0) -- (0,3);
  	\draw (0.7,1) -- (1,0);
  	\draw (1.3,1) -- (1,0);
  	\draw (1.7,1) -- (2,0);
  	\draw (2.3,1) -- (2,0);
  	\draw (3,1) -- (3,0);
  	\draw (0.7,1) -- (1,1.7);
  	\draw (1.3,1) -- (1,1.7);
  	\draw (1,1.7) -- (1,3);
  	\draw (0,2) -- (1,1.7);
  	\draw (0,2) -- (1,2.3);
  	\draw (1.7,1) -- (1.7,2);
  	\draw (1.7,2) -- (2,3);
  	\draw (2.3,2.3) -- (2.3,1);
  	\draw (2.3,2.3) -- (2,3);
  	\draw (1.7,2) -- (1,1.7);
  	\draw (1.7,2) -- (1,2.3);
  	\draw (1.7,2) -- (2.3,2.3);
  	\draw (1.7,2) -- (2.3,1.7);
  	\draw (3.3,2.3) -- (2.3,2.3);
  	\draw (3.3,1.7) -- (2.3,1.7);
  	\draw (3.3,2.3) -- (3.3,1.7);
  	\draw (2.7,2.3) -- (2.7,1.7);
  	\draw (2.7,1.7) -- (3,1);
  	\draw (3.3,1.7) -- (3,1);
  	\draw (2.7,2.3) -- (3,3);
  	\draw (3.3,2.3) -- (3,3);

  	\draw[->,line width=3pt] (0,3) -- (0,2);
  	\draw[->,line width=3pt] (0,2) -- (1,1.7);
  	\draw[->,line width=3pt] (1,1.7) -- (1.3,1);
  	\draw[->,line width=3pt] (1.3,1) -- (1,0);
  	\draw[->,line width=3pt] (1,0) -- (2,0);
	\draw[->,line width=3pt] (2,0) -- (2.3,1);
  	\draw[->,line width=3pt] (2.3,1) -- (3,1);

  	\node at (0,0) {};
  	\node at (1,0) {};
  	\node at (2,0) {};
  	\node at (3,0) {};
  	\node at (0,1) {};
  	\node at (0.7,1) {};
  	\node at (1.3,1) {};
  	\node at (1.7,1) {};
  	\node at (2.3,1) {};
  	\node at (3,1) {};
  	\node at (0,2) {};
  	\node at (1,1.7) {};
	\node at (1,2.3) {};
	\node at (1.7,2) {};
	\node at (2.3,1.7) {};
	\node at (2.3,2.3) {};
	\node at (2.7,1.7) {};
	\node at (2.7,2.3) {};
	\node at (3.3,1.7) {};
	\node at (3.3,2.3) {};
	\node at (0,3) {};
	\node at (1,3) {};
	\node at (2,3) {};
	\node at (3,3) {};
  \end{tikzpicture}
  \caption{The routing during measurement planning prefers regions of low data density. On the right: Corresponding path in the induced abstract graph.}
\end{figure}

\subsubsection*{Measurement and symmetric grid refinement}

The actual measurement is combined with an adaptive refinement of the grid.
The goal of the iteration is to fill the solution map in a way that the remaining small holes can be closed by extrapolation of the surrounding data.
The next step is then:
\begin{enumerate}[I.]
\setcounter{enumi}{2}
\item \emph{Measurement:}
  The actuator settings are varied continuously along the path prescribed by the measurement ramps.
  This results in a linear ordering of the measurements.
\end{enumerate}

For the recording of the $\ell$ measurements per ramp, several additional aspects have to be taken into account, as explained below.
It is important not to store every measurement, since we do not want to store redundant data.
Thus, we restrict our attention to \emph{relevant value variations}, i.e., we only store the measurement $F(u^q)$ at a point $u^q$ if the measurement is sufficiently different.
Consequently, fewer than $s\cdot \ell$ data points may be added to the set $\dmeasurements$.
Here one can choose among various meaningful distance functions.
However, we also want to establish a \emph{minimal measurement frequency}.
That is, if we omitted too many subsequent measurements due to the previous rule, then we store the measurement nonetheless.

Moreover, it may happen that a measured value lies outside the noncritical range $\codomain$, i.e., it violates one or more restrictions.
In order to \emph{exclude critical values}, the entire measurement is disrupted, and we continue from scratch at \eqref{it:random-grid-box} with the last valid measurement.
Special care is needed for observables with a pronounced latency.
The development of these values is extrapolated, and the measurement is rerouted already if the measured values get sufficiently close to the boundary of $\codomain$.

\subsection{Data Cleaning}\label{sec:data-cleaning}
\label{subsec:cleaning}

Assume that the measurement phase of the iteration step is complete, i.e., a new set
\[
\dmeasurements=\{(u^1,y^1),\dots,(u^{M'},y^{M'})\}
\]
is available. In Section \ref{sec:grid-refinement}, we will compute local polynomial fits to the function $F$ using data points in $\dmeasurements$ as interpolation points. At points in which the interpolation is not good enough, the grid is refined further.
A fundamental necessity for the good fit of a polynomial approximation is that the interpolation points from which it is generated resemble a random point cloud. This is not the case for $\dmeasurements$, since all points lie on a piecewise linear path given by the sequence of measurement ramps. We thus need to extract a subset $\dmeasurements_{\text{red}}$ from $\dmeasurements$ that is sufficiently generic.

Furthermore, we want to translate the discrete version \eqref{eq:discrete-opt} of our optimization problem into an integer linear program which picks the right kind of measurements from which we can then obtain our engine maps; cf.\ Section~\ref{sec:discrete-optimization} below.
However, our measurements need to be preprocessed in order to make such an approach feasible by significantly reducing the amount of data without affecting the accuracy.

As it turns out we can achieve both goals with a single method, which we call the \emph{adaptive space compressor}.
The idea is the following.
For the set $\widehat U = \{ u^1,u^2,\dots,u^M \}$, which is the projection of the set $\dmeasurements$ onto the admissible domain $\domain$, we define a corresponding \emph{threshold graph}.
Its nodes are the points $u^q$, and there is an edge between $u^q$ and $u^r$ if the Euclidean distance $\norm{u^q-u^r}_2$ is below a certain threshold.
This threshold depends on the sizes of the two grid boxes which contain $u^q$ and $u^r$.
Then we iteratively remove measurements which have the maximal node degree in the threshold graph until only isolated nodes remain.
This \emph{reduced} set is denoted as $\widehat U_{\text{red}}$, which yields
\[
\dmeasurements_{\text{red}} = \bigl\{ (u,y) \in \dmeasurements \suchthat u \in \widehat U_{\text{red}} \bigr\}.
\]
The set $\dmeasurements_{\text{red}}$ has significantly fewer data points compared to $\dmeasurements$. But due to the criteria for grid refinement, which will be introduced in (V) and (VI) below, the distribution of the remaining data points reflects the structure of the grid $G$, which in turn reflects the behavior of $F$.
This especially means that the data density is higher in regions where the behavior of $F$ is hard to reconstruct by interpolations. Thus, the adaptive space compressor filters out the relevant information gained by the preceding measurements. Moreover, it breaks up the piecewise linear structure of the data point distribution, leading to local (pseudo-)randomness which is required by the polynomial interpolation.
\begin{remark}
  In our experiments we further observed that this approach has the additional benefit of producing solutions which are stable with respect to hysteresis effects.
\end{remark}
Note that this reduction does not affect the set $\dmeasurements$, but we rather explicitly keep $\dmeasurements_{\text{red}}$ as a subset.
In this way, all our measurements are taken into account for computing the weights according to \eqref{eq:wB} and \eqref{eq:wBB} in (I) in the next round of measurements.

\subsection{Grid Refinement}\label{sec:grid-refinement}

The purpose of refining the grid $\grid$ is to accumulate sufficiently many ``meaningful'' data points in the set $\dmeasurements_{\text{red}}$, such that a solution map of $k^2$ data points can be extracted from it which conforms to both emission and drivability constraints.
Recall that in our application we have $m\geq 8$ actuators and $n\geq 14$ measurands. This makes computing on a uniform grid infeasible.
The price to pay is that it is slightly more involved to determine points $u^q$ in the admissible domain $\domain$ such that the data point $(u^q,y^q)$ adds significant information to our data set $\dmeasurements_{\text{red}}$ and is thus stored.

We proceed with the following step:
\begin{enumerate}[I.]
\setcounter{enumi}{3}
\item \emph{Computation of local fit:}
  We employ a Newton interpolation approach, which requires only a partial recomputation of the interpolation polynomial if some interpolation points are exchanged. Each component function $[F]_i: \Rm\rightarrow \R$, where $i\in[n]$, has to be interpolated by a Newton polynomial $L_d[F]_i$, where $d$ is the approximation order.
  Assuming the $(d+1)$-times continuous differentiability of $[F]_i$ at $u$, a $d$-th order polynomial fit requires $N \define {{m+d}\choose{m}}$ interpolation points $\{u^{q_1}, u^{q_2}, \dots, u^{q_N}\}\subseteq \umeasurements$.
  These $N$ points have to satisfy certain spatial relations for the approximation properties of the interpolation polynomials to hold; the interpolation problem is then called \textit{poised}.
  The formulation of the latter conditions is somewhat technical.
  We thus omit the details and refer to the literature instead; cf., e.g., \cite{Schueler2000}.
  Our data cleaning method in Section \ref{sec:data-cleaning} takes care of poisedness.

  If an interpolation polynomial with $N$ interpolation points which have the required geometrical structure fails to give a $d$-th order approximation of $[F]_i$ at $u$, this indicates insufficient smoothness of $[F]_i$ at $u$, which we will use as a criterion for a local grid
  refinement
  We aim at a second order fit, i.e., $d=2$, which requires
  \[
  N\ =\ {{m+2}\choose{m}}\ =\ \frac{m\,(m+1)}{2}
  \]
  interpolation points. Our goal is to compare the measured value $y^q\in\codomain$ at a point $u^q\in\domain$, where $(u^q,y^q)\in\dmeasurements_{\text{red}}$, with a polynomial interpolation of $F$ in $u^q$ using elements of $\dmeasurements_{\text{red}}$ as interpolation points.
  The error of such an approximation is minimized if we use the $N$ closest neighbors $\{u^{q_1}, u^{q_2}, \dots, u^{q_N}\}$ of $u^q$ in $\domain$ as interpolation points, such that \[\{(u^{q_1},y^{q_1}), (u^{q_1},y^{q_2}), \dots, (u^{q_N},y^{q_N})\}\subseteq\dmeasurements_{\text{red}}\,.\]
\end{enumerate}

Concerning the implementation of the polynomial interpolation: As can be seen in the experimental section below, we will usually have less than $50,000$ points in $\dmeasurements_{\text{red}}$. This makes the following brute-force approach feasible.
Let $M_{\text{red}} \define \card{\dmeasurements_{\text{red}}}$ and assume for simplicity that $q=1$.
Now calculate the squared Euclidean distances of $u^2, u^3, \dots, u^{M_{\text{red}}}$ to $u^1$. This has a cost of roughly $3\cdot m\cdot M_{\text{red}}$ elementary arithmetic operations. Store these values in the first row of a $(2\times [M_{\text{red}}-1])$-array and the indices of the corresponding points $u^r$ in the second row ($\approx 2\cdot M_{\text{red}}$ writes). Then determine the $N$ smallest entries of the first row and return the corresponding second row entries.

A simple in-place algorithm to accomplish this is the following: First, determine the $N$-th smallest element of the first row, which costs one sweep of the array; cf \cite[p.183 ff.]{Cormen:2009:IAT:1614191}. Second, sort all columns to the front of the array whose first row entry is smaller than or equal to the $N$-th smallest first row entry. This costs another sweep of the array. Finally, return the first $N$ second row entries.

Repeating this procedure for all $M_{\text{red}}$ points in $\dmeasurements_{\text{red}}$, we arrive at an approximate cost of $(3\cdot m+4)\cdot M_{\text{red}}^2 \in \mathcal O(m\cdot M_{\text{red}}^2)$ elementary operations.
In our setting, for $50,000$ points and $8$ dynamic actuators, this results in a total of about $70$ billion elementary operations, which is a fairly insignificant task for modern computers.
For a general discussion on the nearest neighbor search in high dimensions we refer to the survey \cite{AndoniIndyk:2018}.

The next step is the following:
\begin{enumerate}[I.]
  \setcounter{enumi}{4}
\item \emph{Symmetric grid refinement:}
  To decide whether the grid needs to be refined, we compare the interpolated value, say $\tilde F(u^q)$, with $y^q=F(u^q)$.
  If the deviation exceeds a threshold, then the box $B$ which contains~$u^q$ is split into $2^m$ smaller grid boxes, symmetrically in all coordinate directions.

  The quality of the local fit depends on the differentiability properties of the approximated function. Hence, the symmetric grid refinement increases the measurement density in areas of potential nondifferentiability.
\end{enumerate}
In practice, the static actuators can be ignored for the computation of the local fit and the symmetric grid refinement.

\subsubsection*{Making the Grid Nonuniform}

Recall that ultimately we want to solve the  discrete version \eqref{eq:discrete-opt} of the constrained optimization problem \eqref{eq:continuous-opt}, i.e., we want to minimize the fuel consumption subject to the drivability and emission constraints.
For the iteration step, the drivability constraints are irrelevant.
Therefore we focus on the subset $\emission$ of measurands corresponding to the emissions.
The final output of the calibration will be a set $\SOL$ of $k^2$ data points $d_{ft}$, where $f,t\in[k]$, which cover the $k$-operation field,  while satisfying the emission constraints \eqref{eq:em-constraint-disc}.
In practice, it is a major challenge to find sufficiently many points which satisfy the conditions \eqref{eq:em-constraint-disc} imposed by emission control.
This leads us to a second type of grid refinement.

\begin{enumerate}[I.]
  \setcounter{enumi}{5}
\item \emph{Asymmetric grid refinement:}
  Consider the point $u^M\in\domain$ at which the last measurement was performed.
  Let $f$ and $t$ be such that $(u^M_{\frequency},y^M_{\torque})\in\OP_{ft}$ for $y^M=F(u^M)$.
  If, for any $p\in\emission$, we have
  \[ [F(u^M)]_p\ < \ \min_{(u^q,y^q)\in S_{ft}}([F(u^q)]_p),\]
  i.e., if the emission measurement $[F(u^M)]_p$ is lower than any other value on the corresponding stack $S_{ft}\subset\dmeasurements$, a \emph{cross-measurement} is performed.
  To this end, each actuator is varied individually to determine the direction with the biggest impact on $[F]_p$.
  Afterwards, the grid box $B$ containing $u^M$ is split into two congruent (sub-)boxes
  along the axis corresponding to the actuator whose variation has the biggest impact on $[F]_p$.
\end{enumerate}
Both types of grid refinement do not add data points to $\dmeasurements_{\text{red}}$ directly.
However, the grid refinements increase the probability for picking one of the subboxes in (I).
This way one can hope to find data points near $(u^q,y^q)$ that add relevant information about $F$ and near $(u^M,y^M)$ with better emission values than those currently stored in $\dmeasurements$.

\subsection{Filling the Gaps}\label{sec:hunt-mode}

The calibration algorithm traces the behavior of the map $F$, given by the physical test engine, by constructing a sequence of (one-dimensional) measurement ramps through the domain $\domain$ which is high-dimensional, as is the range of $F$.  Naturally, this approach cannot produce a sufficient coverage of $\domain$.  In particular, often several stacks $S_{ft}$ of the $k$-operation field will remain empty.
If an empty stack $S_{ft}$ lies in the interior of the operation field, we
pick $N=m(m+1)/2$ data points $\{(u^{q_1},y^{q_1}), (u^{q_2},
y^{q_2}),\dots,(u^{q_N},y^{q_N})\}$ from neighboring stacks and close the
hole by a Newton type interpolation on $\{u^{q_1}, u^{q_2},
\dots,u^{q_N}\}$ as done in the local fit step IV.
However, if an empty stack $S_{ft}$ lies on the boundary of the operation field, then the holes are closed by calculating secants of neighboring data points in the nearest stacks and extending them linearly.

\subsection{An Integer Linear Program}\label{sec:discrete-optimization}

In this section we detail how to employ an integer linear program to derive a solution $\SOL$ from the reduced set of measurements $\dmeasurements_{\text{red}}$ that satisfies the objectives listed in Section \ref{sec:discretization}.
For an introduction to the solution of integer linear programs, see \cite{Cook1998}, \cite{Schrijver1998}, and \cite{Sierksma2001}.

Now, for a data point $d^q$, let $S^q$ be the stack that contains it.
As before, let ``$\fuel$'' be the data point index corresponding to the injected fuel quantity.
Then the \emph{prey value} of $d^q$ is defined as
\[
p^q \define \frac{1}{[d^q]_\fuel} \min\bigl\{[d^r]_\fuel \suchthat d^r\in S^q \bigr\}.
\]
For each index $q\in[M]$, where $M \define \card{\dmeasurements_{\text{red}}}$, we introduce a binary decision variable $s^q\in\{0,1\}$, which indicates whether the data point $d^q$ is part of the solution $\SOL$.
Then the objective function of our integer linear program is defined as
\begin{equation}\label{eq:opt}
  \text{minimize}\  \sum_{q=1}^M p^q\, s^q.
\end{equation}

The solution $\SOL$ should contain at most one element from each stack $S_{ft}$.
This condition is reflected in the \emph{stack constraint}
\[
s_{ft}\, +\, \sum_{d^q\in S_{ft}}s^q\ =\ 1 \quad \text{for all } (f,t)\in[k]\times[k],
\]
where $s_{ft}$ is a \emph{stack decision variable}; it is $1$ if it is not possible to choose a data point for stack $S_{ft}$ and $0$ otherwise. The $s_{ft}$ variables will be used in~\eqref{eq:integral-emission-constraint} below to assign penalty values to stacks whose data points are not part of the solution.

To avoid engine damage, the variation speeds of all $m$ actuators are constrained individually by nonnegative constants $\Delta_\actuator$ for $\actuator\in[m]$.
Let $d^q\in S_{ft}$, and $d^r$ be a data point in either neighboring stack $S_{f\pm 1,t}$ or $S_{f,t\pm 1}$.
As our \textit{drivability constraints} we require
\begin{align}\label{eq:drivability-constraints}
s^q\ +\ s^r\ \leq\  1\qquad\text{if}\quad \bigl\lvert [d^q]_\actuator-[d^r]_\actuator \bigr\rvert \geq \Delta_\actuator\quad \text{ for any } \actuator\in[m].
\end{align}
Note that this is a secant constraint, since the data points compared are contained in neighboring stacks.
The emission condition \eqref{eq:em-constraint-disc} describes the upper bound of several emission test cycles, e.g., maximal \NX production.
Additionally, the current output of a pollutant $p\in\emission$ is restricted by the boundaries of the respective interval of noncritical values
\[
J_p = [\underline{e}_p,\, \overline{e}_p].
\]
In the following, the weights $\omega_{ft}$ represent the mean resistance time of the of the rectangles $\OP_{ft}$ in the given test cycle.
For simplicity of notation, the mean resistance time on the rectangle corresponding to the stack $S^q$ containing $d^q$ is denoted by $\omega^q$.
Then the \emph{emission constraint} can be formulated as follows:
\begin{align}\label{eq:integral-emission-constraint}
  \sum_{q=1}^{M} \omega^q\, [d^q]_p\, s_q + \sum_{(f,t)\in [k]\times [k]} \omega_{ft}\, \overline{e}_p\, s_{ft}\ \leq\ \emconstraint_p\qquad\text{for all } p\in\emission\,.
\end{align}
The second summand of the left hand side of the inequality serves as a penalty term to compensate for stacks that contribute no data point to the solution $\SOL$.

\section{AVL Engine Model}\label{sec:avl}

For our experiments in Section \ref{sec:experiments} we replace the test-bench with an engine model of a diesel engine with turbo charger, pilot injection and variable turbine geometry.
The latter has been developed in cooperation with AVL GmbH and is based on measurements on a compression ignition/diesel engine.
Below we will give a brief overview of the effects of the $8$ actuators and $18$ measurands that are simulated, thus indicating the scope of the simulation.
For a detailed description of the AVL model's derivation, see \cite{Burggraf} and~\cite{AVL}.

\subsection{Actuators of the AVL Engine Model}\label{sec:actuators}

\paragraph{\it Revolution frequency of the crank shaft (\RF)}
The crankshaft converts the reciprocating motion of the cylinders into a rotational motion.
In modern four-stroke engines every cylinder fires once for every two revolutions of the crankshaft.
The revolution frequency is a controlled actuator, as discussed in Section \ref{sec:math}.
It stands out among the other actuators since it provides one coordinate axis of the operation field.

\paragraph{\it Injected fuel quantity (\IF)} In contrast to a spark-ignited engine, the injected amount of fuel \textit{Q} is the most important actuator. More precisely, the injection process is crucial in the application process of diesel engines. The injection process is given by several pre/pilot-injections, a main injection and post-injections. Typically, the engine torque is mainly determined by the main-injection \textit{Q}. This actuator defines the total amount of fuel per cycle. In this simple model, the injected fuel volume is divided into one pilot and the main injection.

\paragraph{\it Pressure in the common rail system (\RP)} In contrast to solenoid-controlled unit injector elements, the pressure is generated by a central fuel and high pressure pump. The fuel injectors are opened and closed by piezo elements.

\paragraph{\it Air filling (\AF)} Similar to a spark ignited engine, an air valve controls the amount of air which contributes to the combustion process. In this model the amount of air is given directly in mass per piston stroke.

\paragraph{\it Turbine geometry (\TG)} Modern turbochargers do not have a static turbine geometry. Variable-turbine-geometry turbochargers are able to tune the angle of the turbine blades in order to increase the amount of boost. Alternative setups are given by static turbines with waste gates. Waste gates are applied to reduce the amount of exhaust gas that accelerates the turbine, so the amount of boost can be controlled. In our model the geometry is given as a value between $30$ and $85$.

\paragraph{\it Main timing (\MT)} The main timing is comparable to the spark timing of Otto-engines. It defines the start timing of chemical reactions in the crank\-shaft angle of the main injection. Similar to the Otto engine, the pressure rise is delayed by the ignition delay.

\paragraph{\it Pilot injection (\PI) and pilot timing (\PT)} Pilot injection works in tandem with pilot timing to achieve a complete burning of the fuel, which in turn also drastically reduces the emission of \NX gases. The pilot injection increases the temperature of the combustion chamber, thus when the main injection occurs the fuel is sent into a chamber which already is at a higher temperature than its autoignition point. This especially facilitates the fuel burning at lower speeds.

\renewcommand{\arraystretch}{1.3}
\begin{table}
  \caption{\label{tab:actuators} Actuator-intervals during various calibration runs. LT/HT = low/high torque, FR = free variation of torque; cf.\ Sections \ref{sec:basic-calibration} and \ref{sec:full-calibration}. }
	\small
	\centering
	\begin{tabular*}{\textwidth}{@{}l@{\extracolsep{\fill}}lllll@{}}\toprule
		Act. & Unit & Basic, LT & Basic, HT & Basic, FR & Full Calibration \\ \midrule
		 \RF & $\frac {1}{\text{min}}$ & 1000--2600 & 1000--2600 & 1000--2600 & 1000--2600 \\
		 \IF & $\frac{\text{mm}^3}{\text{cycle}}$ & 6--10 & 50--60 & 6--60 & 6--60 \\
		 \RP & $\operatorname{hPa}$ & 295677 / 405677 & 295677 & 295677 & 295677--1126537 \\
		 \AF & $\frac{\text{mg}}{\text{stroke}}$ & 300 & 300 & 300 & 275--991 \\
		 \TG & $\operatorname{int}$ & 30 & 30 & 30 & 30--85 \\
		 \MT & $ ^\circ \operatorname{CA}$ & 0 / 10 & 0 & 0 & 0--10 \\
		 \PI & $\frac{\text{mm}^3}{\text{cycle}}$ & 1 & 1 & 1 & 1 \\
		 \PT & $\mu s$ & 1540 & 1540 & 1540 & 1540--2565	\\\bottomrule
	\end{tabular*}
\end{table}

\subsection{Measurands of the AVL Engine Model}

\paragraph{\it Torque} The produced torque of the engine. The revolution frequency and the engine torque define the power level of the ICE.

\paragraph{\it Fuel mass flow} The amount of fuel which is delivered to the cylinder per hour.

\paragraph{\it Carbon monoxide, hydrocarbons, nitrogen oxides, soot}  The momentary and integral exhaust emission of pollutants, given in parts per million and emission per hour.

\paragraph{\it Indicated mean pressure} The average pressure over a cycle in the cylinder.

\paragraph{\it Lambda} Lambda displays the chemical partitioning of the fuel. Spark-ignited engines run on lambda around $1$. The combustion limits are 0.6 and 1.6. The lambda value for compression ignition engines is much higher, up to 20.

\paragraph{\it Manifold pressure} The manifold pressure measures the absolute pressure in front of the intake channel. The pressure depends on the absolute environmental pressure and the boost level of the turbocharging system.

\paragraph{\it Boost pressure} The boost pressure is created by the turbocharging unit. It depends on the turbine geometry setting and the current combustion behavior.

\paragraph{\it Maximal cylinder pressure} The maximal point of the cylinder pressure sequence. Every engine has a specified maximal cylinder pressure, in order to avoid damaging of the devices. Typically, the maximal pressure is about $160$ bar.

\paragraph{\it Manifold temperature} The manifold temperature measures the temperature in the intake channel.

\paragraph{\it Critical temperature} The critical temperature is defined as the temperature of the burn zone when the exhaust valves opens. In case of bad timings, the fuel has not been consumed completely, which results in the release of flames to the exhaust manifold, catalysts and turbocharger. It indicates damages to sensitive parts of the engine setup.

\paragraph{\it Specific fuel consumption} The current power level of the engine over the current fuel consumption.

\section{Optimizing the AVL Engine Model}\label{sec:experiments}

We subdivided the engine optimization into two phases. During the first phase, which we dub basic calibration, the emission constraint~\eqref{eq:em-constraint} is ignored. The purpose of this omission is to obtain a coarse picture of the engine behavior and derive a preliminary solution that conforms to the drivability constraint \eqref{eq:drivability-constraints}. In the second phase, which we dub full calibration, a solution conforming to all constraints is obtained by building on and refining the preliminary solution of the first phase.

All measurement times are given as real-world measurement times as opposed to runtimes of the simulation which are highly hardware specific.

\subsection{Basic Calibration}\label{sec:basic-calibration} During basic calibration, only two actuators are dynamic, \IF and \RF. The other six actuators are static throughout.
The phase is subdivided into three runs of several hours. By a ``run'' we mean the following:
The intervals of the dynamic actuator \IF and the values of the static actuators are reset. Then the algorithm steps in Sections~\ref{sec:iteration-step}--\ref{sec:grid-refinement} are repeated for a prescribed amount of time.
The revolution frequency is set to vary over its full range of 1000--2600 revolutions per minute throughout all three runs.
For a full account of the applied settings, see Figure \ref{tab:actuators}.

A good foundation for an optimization is a well-defined low and high torque boundary region. In a compression ignition engine the generated torque is roughly proportional to the injected fuel. In the first run we measure the low torque region by limiting the \IF-interval to 6--10 $\text{mm}^3/\text{cycle}$. The static actuators are set to values that support the creation of low torque operation points.
The measurement time of this run is $6$ hours. The low torque region is measured for a second time span of $6$ hours with slightly modified settings, i.e., \MT $= 10.0$ and \RP $= 405677 \operatorname{hPa}$.

Afterwards, the \IF-interval is reset to 50--60 $\text{mm}^3/\text{cycle}$, so that the high torque regions can be measured.
The static actuators are set to values that support high torque operation points.
Measurement time is again $6$ hours.

Finally, the lower and upper bound of \IF are removed, in order to get a picture of the remaining region of the solution map. Therefore, the static actuators are set to midrange values. This region is measured for $6$ hours, too. After $24$ hours, we get a coarse picture of the engine behavior.

\begin{figure}
	\includegraphics[width=\textwidth]{scattering.pdf}
	\caption{Operation field saturation after basic calibration, before closing of gaps.}
\end{figure}

Using the data as input for an integer linear program as described in Section \ref{sec:discrete-optimization} yields a solution map that still has holes. These are caused by the drivability constraint.
These are closed by interpolation/extrapolation as described in Section \ref{sec:hunt-mode}.
For the comparison of a solution that respects the drivability constraint with one that does not, see Figure \ref{fig:af-comp}.

\begin{figure}[htbp]
	\begin{center}
		\begin{minipage}{0.9\textwidth}
			\centering
			\includegraphics[width=1.0\textwidth]{stat_sol_RANDOM_air.pdf}

		\end{minipage}
	\end{center}
	\vfill
	\begin{center}
		\begin{minipage}{0.9\textwidth}
			\centering
			\includegraphics[width=1.0\textwidth]{Air_Filling_RANDOM.pdf}
		\end{minipage}
	\end{center}
	\caption{Top: Solution for actuator \AF without drivability constraint. Bottom: \AF settings of \euro 5 solution respecting the drivability constraint.}
	\label{fig:af-comp}
\end{figure}

In the current example run it took an additional $45.8$ hours to close all measurement holes. The average fuel consumption for $100$ kilometers of the so-derived preliminary solution is $4.65$ liters.
Since this optimization does not take exhaust emissions into account, the cycle integrals of \CO, \HC and \NX are rather high: 3.68, 1.72 and 5.26, respectively, per \nedc. The solution map obtained in the first calibration phase serves as the basis of the full calibration.

\subsection{Full Calibration}\label{sec:full-calibration}

During the full calibration phase the only static actuator is pilot injection. All others range over their whole respective domains. There is no predefined time limit for the iteration. The algorithm steps in Sections~\ref{sec:iteration-step}--\ref{sec:grid-refinement} and the solution of the integer linear program in Section~\ref{sec:discrete-optimization} are repeated until a solution map conforming to the given constraints is derived.

The optimization is performed, first, for the \nedc, then for the \rand cycle. An \nedc solution that conforms to the \euro 4 norm is obtained after $68.51$ hours, one that conforms to \euro 5 after $88.41$ hours. From the data set of the \nedc \euro 5 solution, a \euro 4 solution for the \rand cycle can be derived.
A \euro 5 solution is obtained after $107.55$ hours of measurement. (For the emission constraints corresponding to the different \euro norms, see Table \ref{tab:euro-norms}.)
To lend some context to these numbers, we performed the same calibrations
using uniform grids. To derive a solution of equivalent quality, i.e., a solution conforming to all constraints given by the different \euro norms, the grid approach consistently required 15 to 20 times as many measurements.

Since the AVL model only has a measurand for particle mass, but not for the number of particles, the determination of a \euro 6 solution lies outside of the model's scope. However, the step from \euro 5 to \euro 6 merely adds an item to the list of emission constraints. This poses no fundamental challenge to our method which is scalable with respect to the number of pollutant limits. (Naturally, adding further constraints will increase the measurement time though.)

\subsubsection*{Description of Figures}

Figures \ref{fig:nox-nedc} and \ref{fig:nox-random} illustrate a typical feature of engine calibrations fitted to specific driving cycles. These essentially subdivide the solution map into two parts. One covered by the cycle where emissions are optimized, and one where they are not. This leads to emission values being 20--100 times higher on the part of the solution map that is not covered by the driving cycle. The \nedc covers only a fraction of the
operation field, while the \rand cycle covers more than half. This results in
lower \NX output on the better part of the operation field for the \rand calibration; cf Figure \ref{fig:diff-nox}.
Conversely, in Figure \ref{fig:diff-fl} one can observe that the calibration solution for the \rand cycle displays 10--15\,\% higher values for specific fuel consumption on most points of the operation field than the corresponding operation points of the \nedc.

One can, of course, enforce the selection of data points with lower emission values on the part of the operation field not covered by a driving cycle as well, but at the expense of a higher fuel consumption.

\begin{figure}[htbp]
	\begin{center}
		\begin{minipage}{0.9\textwidth}
			\centering
			\includegraphics[width=1.0\textwidth]{Weight_NEDC.pdf}
		\end{minipage}
	\end{center}
	\vfill
	\begin{center}
		\begin{minipage}{0.9\textwidth}
			\centering
			\includegraphics[width=1.0\textwidth]{NOx_NEDC.pdf}

		\end{minipage}
	\end{center}
	\vfill
	\begin{center}
		\begin{minipage}{0.9\textwidth}
			\centering
			\includegraphics[width=1.0\textwidth]{NOx_NEDC_log.pdf}

		\end{minipage}
	\end{center}
	\caption{Top: Weighting table of the \nedc. Middle: Resulting \NX emissions (in $g/h$) after calibration for the \nedc. Bottom: \NX emissions with logarithmic scale.}
	\label{fig:nox-nedc}
\end{figure}

\begin{figure}[htbp]
	\begin{center}
		\begin{minipage}{0.9\textwidth}
			\centering
			\includegraphics[width=1.0\textwidth]{Weight_RANDOM.pdf}
		\end{minipage}
	\end{center}
	\vfill
	\begin{center}
		\begin{minipage}{0.9\textwidth}
			\centering
			\includegraphics[width=1.0\textwidth]{NOx_RANDOM.pdf}

		\end{minipage}
	\end{center}
	\vfill
	\begin{center}
		\begin{minipage}{0.9\textwidth}
			\centering
			\includegraphics[width=1.0\textwidth]{NOx_RANDOM_log.pdf}

		\end{minipage}
	\end{center}
	\caption{Top: Weighting table of the \rand cycle. Middle: Resulting \NX emissions (in $g/h$) after calibration for the \rand cycle. Bottom: \NX emissions on a logarithmic scale.}
	\label{fig:nox-random}
\end{figure}

\begin{figure}[htbp]
	\begin{center}
		\begin{minipage}{0.9\textwidth}
			\centering
			\includegraphics[width=1.0\textwidth]{diff_NEDC_RANDOM_NOx.pdf}
		\end{minipage}
	\end{center}
	\caption{\footnotesize
		Operation points where calibration for \rand cycle yields higher \NX emission than in \nedc. Colors represent the ratio of the emission values. }
	\label{fig:diff-nox}
\end{figure}

\begin{figure}[htbp]

	\begin{center}
		\begin{minipage}{0.9\textwidth}
			\centering
			\includegraphics[width=1.0\textwidth]{Fuel_NEDC.pdf}

		\end{minipage}
	\end{center}
	\vfill
	\begin{center}
		\begin{minipage}{0.9\textwidth}
			\centering
			\includegraphics[width=1.0\textwidth]{Fuel_RANDOM.pdf}

		\end{minipage}
	\end{center}
\vfill

	\begin{center}
		\begin{minipage}{0.9\textwidth}
			\centering
			\includegraphics[width=1.0\textwidth]{diff_NEDC_RANDOM_fl.pdf}

		\end{minipage}
	\end{center}
	\caption{Actuator \IF after calibration for the \nedc (top) and \rand cycle (middle). Bottom: Operation points with higher values in calibration for \rand cycle than for \nedc. Colors represent ratio of consumption values.}
	\label{fig:diff-fl}
\end{figure}

\section{Conclusion}

In this article, we have described a semi-automatic approach to
calibrate and optimize internal combustion engines with several
actuators and sensors. The side constraints are limits given by
safety or technical requirements, bounding the variation speed of actuators and
ensuring emission bounds on given driving cycles. Our method automatically
performs refinements of measurements, thus focusing the effort of the measurements around
regions of nonsmooth behavior of the engine. That is, it automatically identifies neighborhoods
of anomalous engine behavior and maps them in appropriate detail.
For the so-obtained data an optimal calibration solution is computed.

The output of the algorithm is a solution map which consists of actual measurements
and thus reflects the exact behavior of the engine for the given settings, as opposed
to indirectly derived actuator setting obtained via modeling or interpolation.
This results especially in improved values for pollutant emission and fuel consumption near nonsmooth regions of the admissible domain.

In our experiments, we demonstrated the practicability of the adaptive meshing methodology,
showing a significant speed-up in the measurement time in comparison to uniform grids.
Moreover, the resulting solutions respect the emission constraints of \euro 4 and 5 norms.
We would like to stress that in our method it is easy to take into account further emission constraints such as, e.g., the number of emitted particles for the \euro 6 norms.

The next interesting step will be to test the method on an actual combustion engine. It is
our expectation that the experimental findings of this work will transfer well to the
real-life setting which is the engine test bench. That is, we expect the concept of adaptive
grid refinement to considerably improve on uniform grids in terms of the number of measurement
required to obtain a sufficiently meaningful image of the engine behavior to derive calibration
settings respecting a given set of constraints from the obtained data.

\bibliographystyle{alpha}
\bibliography{engine}

\end{document}

config: arxiv


TITLE: Is there a quick way of determinging the length of a piece of felt that is needed to cover a space having equal size hills and valleys?
QUESTION [1 upvotes]: In the making of a pen case for expensive fountain pens, I affix equally-spaced wooden half-rounds to the bottom inside surface of, say, a cigar box and cover the surface with felt. The half-rounds and felt prevent the pens from rubbing up against each other. When completed, the case is a good way to display the pens when the hinged top is open, as well as to keep dust out when the hinged top is closed. Jewelers might do something similar in a section of a glass-enclosed display case to display longer items such as bracelets, necklaces, and more.
One obvious way of measuring a piece of felt that will cover completely the hills and valleys of the surface, when cut, is to take a longer-than-needed piece of felt of the correct width, cover the "lands and grooves" with it, and mark the place on the felt where it completely covers the surface. Then cut it and paste it (permanently) to the bottom of the box and to the half-rounds.
As an alternative, however, is there a way to measure the length of felt needed by using a formula for any size of box?
For example, if the bottom of the box is five inches from front to back; if the half rounds have a diameter of 5/8 of an inch; and if the spaces between the half-rounds are the same (say, half an inch, plus or minus); what would the formula be to measure and cut the piece of felt?
The only tag I could think of is "arc-length," as I am not a mathematician, as you can probably tell. By the way, I may in the near future add to this question a two-dimensional drawing of what I'm describing, for your ease in answering my question.
Thanks to an architect friend of mine:

Is there a formula in there somewhere?

REPLY [1 votes]: If I understand the question correctly, it's pretty simple.  If you have $n$ segments where each segment consists of a flat area of width $w$ followed by a semicircle of diameter $d$, then the total arc length is
$$s = n(w + \pi d/2)$$
In the drawing we have
$$\begin{align}
n &= 5\\
w &= 1/2 \\
d &= 5/8
\end{align}$$
so
$$s = 5 \times (0.50 + 3.14 \times 0.625 / 2) = 5 \times 1.481 = 7.41$$


text:
\section{Proof of Theorem \ref{thm: other direction (intro)}}
\label{sec: the other direction}

\begin{definition}
	Let $G$ be a group, and $S\subseteq G$ a subset. A subgroup $H \le G$ is \emph{$S$-conjugacy complete} if it intersects the conjugacy classes of all elements of $S$.

	When $S=G$ we say that $H$ is \emph{conjugacy complete}.
\end{definition}

It was observed in \cite{kantor2011invariable} that the following are equivalent definitions of IG.
\begin{lemma}\label{lem: equivalent definitions of IG}
	Let $G$ be a group, and $S\subseteq G$ a subset.
	The following are equivalent:
	\begin{enumerate}
	    \item $S$ invariably generates $G$ \label{equivalent definitions - IG}
	    \item $G$ does not contain a proper $S$-conjugacy complete subgroup.\label{equivalent definitions - conjugacy complete subgroup}
	    \item Every non-trivial transitive action $G\actson X$ has an element $s\in S$ without fixed points.\label{equivalent definitions - transitive action}
	\end{enumerate}
\end{lemma}

Wiegold \cite{wiegold1976transitive}
proved that the class of IG groups is closed under extensions, in fact the following slightly stronger result holds.

\begin{proposition}\label{prop: other direction}
Let $G$ be a group,

$N\le H \le G$ be subgroups and $N\normalin G$.
Let $S\subseteq H$ and $S'\subseteq G$. If $H$ is invariably generated by $S$ and $G/N$ is invariably generated by $S'$ then $G$ is invariably generated by $S\cup S'$.
\end{proposition}

\begin{proof}
    Let $G\actson X$ be a transitive action on a set with $|X|\ge 2$. We want to find an element of $S\cup S'$ which acts without fixed points on $X$.

    Since $N$ is normal, we know that $G/N\actson X/N$. If $|X/N|\ge 2$ then since $G/N$ is invariably generated by $S'$, there exists an element $s' \in S'$ that acts without fixed points on $X/N$ and hence also on $X$. If $|X/N|=1$, then $N$, and hence $H$, act transitively on $X$ and since $H$ is invariably generated by $S$, there is an element $s\in S$ which acts without fixed points on $X$.
\end{proof}

In particular, we can deduce Theorem \ref{thm: other direction (intro)}.
\begin{proof}[Proof of Theorem \ref{thm: other direction (intro)}]
If $H$ is a finite index IG (resp. FIG) subgroup of $G$, then $N=\Core_G(H)=\bigcap_{g\in G} g\ii Hg$ is of finite index in $G$. Since every finite group is FIG, we get that $N\le H\le G$ satisfy the assumptions of Proposition \ref{prop: other direction} which implies that $G$ is IG (resp. FIG).
\end{proof}


TITLE: Series termination
QUESTION [0 upvotes]: The successive terms in a power series are given by the recurrence relation $$\frac{a_{n+1}}{a_{n}}= \frac{n(n-1)+\lambda}{9(n+1)(n+2)}$$ where $\lambda=\text{const.}\in\mathbb{R}$. So our power series is \begin{align}\sum^{\infty}_{n=0}a_{n}x^{n} &= a_0 + a_1x + a_2x^2+a_3x^3+\cdots\\&=a_{0}\left(1 + \frac{\lambda}{18}x^{2} + \frac{\lambda(\lambda+2)}{1944}x^{4} + \frac{\lambda(\lambda+2)(\lambda+12)}{524880}x^{6} +\cdots\right) \\ &+ a_{1} x \left(1 + \frac{\lambda}{54}x^{2} + \frac{\lambda(\lambda+6)}{9760}x^{4} \frac{\lambda(\lambda+6)(\lambda+20)}{3674160}x^{6} + \cdots\right) \\ &= a_{0}\,f{(x)} + a_{1}x\,g{(x)}\end{align} Where \begin{align}f(x) &= 1 + \frac{\lambda}{18}x^{2} + \frac{\lambda(\lambda+2)}{1944}x^{4} + \frac{\lambda(\lambda+2)(\lambda+12)}{524880}x^{6} +\cdots \\ g(x) &= 1 + \frac{\lambda}{54}x^{2} + \frac{\lambda(\lambda+6)}{9760}x^{4} \frac{\lambda(\lambda+6)(\lambda+20)}{3674160}x^{6} + \cdots\end{align} I must find all the values $\lambda$ for which either $f(x)$ and $g(x)$ terminate (ie. form a polynomial). Although i'm not really sure what to do. If either one of these terminate then the general term $a_{n} = 0$ for some $n\ge0$. This happens if and only if $\lambda = n(1-n) \quad \forall n\ge0$. But this is using the recurrence relation for both of them together, how would I get a condition for them separately?

REPLY [1 votes]: If lambda equals n(1-n) with n even, then your power series for f(x) terminates. You then get a polynomial for your original power series by choosing a_1 to be 0. (By "original power series" I mean your power series for a_0*f + a_1*x*g).
On the other hand, if lambda equals n(1-n) with n odd, then your power series for g(x) terminates. You then get a polynomial for your original power series by choosing a_0 to be 0.
Thus, assuming neither a_0 nor a_1 is equal to zero, there is no single value of lambda that will yield a polynomial for your original power series!
Lastly, Im not sure where you got this question from, but just so that you know its physical significance: these sort of recurrence relations are very common in solving the Schrodinger equation in quantum theory.
Hope this helps, and please feel free to ask any more questions if you need!