Jong Ho Jhee

Myung Jun Kim

Myeonggeon Park

Jeongheun Yeon

Hyunjung Shin
email: shin@ajou.ac.kr

Fast Prediction for Criminal Suspects through Neighbor Mutual Information-Based Latent Network

One of the interesting characteristics of crime data is that criminal cases are often interrelated. Criminal acts may be similar, and similar incidents may occur consecutively by the same ofender or by the same criminal group. Among many machine learning algorithms, network-based approaches are well-suited to refect these associative characteristics. Applying machine learning to criminal networks composed of cases and their associates can predict potential suspects. Tis narrows the scope of an investigation, saving time and cost. However, inference from criminal networks is not straightforward as it requires being able to process complex information entangled with case-to-case, person-to-person, and case-to-person connections. Besides, being useful at a crime scene requires urgency. However, predictions from network-based machine learning algorithms are generally slow when the data is large and complex in structure. Tese limitations are an immediate barrier to any practical use of the criminal network geared by machine learning. In this study, we propose a criminal network-based suspect prediction framework. Te network we designed has a unique structure, such as a sandwich panel, in which one side is a network of crime cases and the other side is a network of people such as victims, criminals, and witnesses. Te two networks are connected by relationships between the case and the persons involved in the case. Te proposed method is then further developed into a fast inference algorithm for large-scale criminal networks. Experiments on benchmark data showed that the fast inference algorithm signifcantly reduced execution time while still being competitive in performance comparisons of the original algorithm and other existing approaches. Based on actual crime data provided by the Korean National Police, several examples of how the proposed method is applied are shown.

Introduction

As the number of criminal cases continues to rise, there has been a notable increase in the utilization of machine learning (ML) for criminal investigation support. Various areas, including crime pattern analysis [START_REF] Ullah | An intelligent system for complex violence pattern analysis and detection[END_REF][START_REF] Sardana | Supervised and unsupervised machine learning methodologies for crime pattern analysis[END_REF][START_REF] Yan | Research on prediction of multi-class theft crimes by an optimized decomposition and fusion method based on XGBoost[END_REF], fraud detection, trafc violation monitoring, sexual assault investigations, and cybercrime analysis [START_REF] Prabakaran | Survey of analysis of crime detection techniques using data mining and machine learning[END_REF][START_REF] Mcclendon | Using machine learning algorithms to analyze crime data[END_REF][START_REF] Vijayarani | A comprehensive analysis of crime analysis using data mining techniques[END_REF], have seen the application of ML techniques. While challenges related to data confdentiality still exist, the potential of machine learning in aiding case investigations is widely recognized. Te valuable insights derived from accumulated criminal cases play a crucial role in providing clues and assisting law enforcement agencies in their investigative eforts.

In the meantime, one of the characteristics of crime data is that criminal cases are related. Criminal behaviors may be similar, and similar incidents may occur consecutively by the same ofender or by the same criminal group. To study crime cases, social network approaches have emerged that refect the associations between people, and these network-based (or graph-based) methods are expected to become one of the standard tools in criminology research [START_REF] Sparrow | Te application of network analysis to criminal intelligence: an assessment of the prospects[END_REF][START_REF] Bright | Using social network analysis to study crime: navigating the challenges of criminal justice records[END_REF][START_REF] Faust | Social networks and crime: pitfalls and promises for advancing the feld[END_REF]. Tere are early works applied to criminal investigations using networkbased ML algorithms. Weber et al. [START_REF] Weber | Scalable graph learning for anti-money laundering: a frst look[END_REF] applied a graph convolutional network algorithm to fnancial data for antimoney laundering forensic analysis. Te task was to predict the suspiciousness of a given target, and Das et al. [START_REF] Das | A graph based clustering approach for relation extraction from crime data[END_REF] used graph-based clustering to extract relational information from crime data of India. Te named entities were extracted from the text corpus of crimes and converted into vectors using the word2vec algorithm [START_REF] Faust | Social networks and crime: pitfalls and promises for advancing the feld[END_REF]. Te network was then constructed by measuring the similarity between these vectors. Clusters found in the network represent incidents or ofenders with similar patterns. Meanwhile, social networks of online auction users are used for fraud detection. By labeling known scammers and legitimate users, a label propagation algorithm is used to predict potential scammers [START_REF] Bangcharoensap | Two step graph-based semi-supervised learning for online auction fraud detection[END_REF]. Also, in [START_REF] Alvari | Semi-supervised learning for detecting human trafcking[END_REF], an online advertising network was constructed and analyzed using Laplacian SVM to detect human trafcking advertisements. Khan et al. [START_REF] Khan | Predicting and preventing crime: a crime prediction model using san francisco crime data by classifcation techniques[END_REF] proposed a crime prediction model by comparing three known algorithms, Naïve Bayes, random forest, and gradient boosting decision tree, and classifed the top ten crimes from the San Francisco crime data. A combined framework of graph representation learning and machine learning methods is introduced to predict the amount of money exchanged among criminal agents and to recover the missing criminal partnerships [START_REF] Lopes | Machine learning partners in criminal networks[END_REF]. Meanwhile, one of the important parts in criminal network analysis is the link prediction problem. In many cases, there is a possibility of acquiring missing or incomplete information by the crime investigation, and one might want to recover missing links or connections among individuals or resources in the information. In that sense, Berlusconi et al. [START_REF] Berlusconi | Link prediction in criminal networks: a tool for criminal intelligence analysis[END_REF] proposed a method to identify missing links in a criminal network by classifying links based on the topological analysis and applied it to the Italian criminal case data against a mafa group. Also, to make the link prediction robust to varying relations in a criminal network, Calderoni et al. [START_REF] Calderoni | Robust link prediction in criminal networks: a case study of the Sicilian Mafa[END_REF] applied various link prediction algorithms and observed the algorithm that leverages the full graph topology.

Of the many network-based ML algorithms, the graphbased semisupervised learning (GSSL) algorithm is one of the most popular because it is easy to use, can handle situations where data have few labels, and its inference is intuitive along the network structure [START_REF] Zhou | Learning with local and global consistency[END_REF][START_REF] Zhou | A regularization framework for learning from graph data[END_REF][START_REF] Hofmann | Consistency of semi-supervised learning algorithms on graphs: probit and one-hot methods[END_REF][START_REF] Belkin | Manifold regularization: a geometric framework for learning from labeled and unlabeled examples[END_REF][START_REF] Xu | Discriminative semi-supervised feature selection via manifold regularization[END_REF][START_REF] Zhou | Learning from labeled and unlabeled data on a directed graph[END_REF][START_REF] Zhu | Semi-supervised learning using Gaussian felds and harmonic functions[END_REF]. Terefore, the scope of application is wide where relational information is important, such as fnding key genes using disease and gene networks [START_REF] Nguyen | Detecting disease genes based on semi-supervised learning and protein-protein interaction networks[END_REF], predicting protein functions using multiple biological networks [START_REF] Shin | Protein functional class prediction with a combined graph[END_REF], and classifying historical fgures to political parties using many relationships such as blood ties, academic ties, and geographic proximity [START_REF] Lee | Historical inference based on semi-supervised learning[END_REF][START_REF] Kim | Inference on historical factions based on multi-layered network of historical fgures[END_REF]. In the domain of social networks, GSSL is used to create relevant links to the concept referred in Wikipedia for all tweet mentions [START_REF] Huang | Collective tweet wikifcation based on semi-supervised graph regularization[END_REF], and in other applications, it is used to detect fake users from a large volume of Twitter networks [START_REF] Balaanand | An enhanced graphbased semi-supervised learning algorithm to detect fake users on Twitter[END_REF]. In the computer vision domain, GSSL is employed for hyperspectral image classifcation based on image networks [START_REF] Shao | Spatial and class structure regularized sparse representation graph for semisupervised hyperspectral image classifcation[END_REF], and in the natural language processing domain, part-of-speech tagging was performed by applying GSSL on a random feld network [START_REF] Subramanya | Efcient graphbased semi-supervised learning of structured tagging models[END_REF].

In this study, we propose a framework for predicting suspect candidates by applying GSSL to criminal networks. Te network is designed to be layered, like a sandwich panel, with a network of crime cases on one side and a network of people, such as victims, ofenders, and witnesses, on the other side. Nodes or entities in each network are connected through similarities. It is also connected to nodes belonging to the other network refecting the relationships between the case and the person who is involved in the case. Meanwhile, applying GSSL to the crime scene requires agility to achieve immediate results. However, when the data size is large, the time complexity of network-based algorithms increases exponentially according to the size of the network, so the inference speed by GSSL is inevitably slowed down. Given that crime scenes always demand a sense of urgency, the slow inference is fatal indeed to solving cases. Terefore, we propose a fast GSSL algorithm for large-scale criminal networks to mitigate the limitation. Te idea is to insert a latent network of cluster centroids and link cases or persons to the corresponding centroids. Exhaustive searches are avoided by inserting a network between the network of cases and the network of persons. Te scope of the search is reduced to only a small set of nodes (cases or people) belonging to the same cluster in the latent network. On the other hand, high-dimensional data are newly represented as low-dimensional vectors by using the clusters in the latent network. Te proposed method is called neighbor mutual information semisupervised learning (MISSL) because it uses the mutual information between the clusters and the neighbors of the nodes. MISSL is robust to nonspherical clusters of various sizes and shapes. Te number of clusters does not increase linearly with the number of cases or people, which is advantageous in terms of memory efciency, especially for networks that require frequent updates and are constantly growing.

Te remaining sections of the paper are organized as follows. Section 2 describes the method of building the criminal network and explains the prediction procedure using GSSL. Section 3 describes how to insert the latent network and details the fast prediction algorithm, MISSL. Section 4 demonstrates the comparative experiments on benchmarking datasets, and Section 5 presents a practical application of criminal network analysis. Finally, conclusions are drawn in Section 6.

Criminal Network and Suspect Scoring

2.1. Network Construction. Criminal acts often exhibit similarities, and it is not uncommon for similar incidents to occur consecutively, either by the same ofender or by a particular criminal group. In addition, in some cases, the patterns of crimes may bear resemblance, even if the offenders involved are diferent, resembling a copycat scenario seen in serial murder cases. Considering these factors, a network serves as an efective tool for illustrating connections between crimes and people. Te criminal network we have developed captures associations between cases, individuals, and interactions between cases and individuals. Tis network follows a two-layered structure, with one layer representing the network of crime cases and the other layer representing the network of people. Tese two layers are interconnected through relationships that link specifc cases with the corresponding individuals involved.

Let the criminal network be denoted as G � (G p × G c ) where G p stands for the people network and G c stands for the case network. Each network is represented as

G � (V, W) where V � x i 􏼈 􏼉 n i�1 , x i ∈ R d is a set of nodes, and W � w ij 􏽮 􏽯 n i,j�1
is a set of weighted edges. Te weight w ij is determined by the similarity between nodes x i and x j . Generally, the Gaussian kernels are widely used for similarity calculation and are represented as follows:

w ij � exp - x i -x j � � � � � � � � � � 2 2σ 2 ⎛ ⎜ ⎜ ⎝ ⎞ ⎟ ⎟ ⎠ , ( 1
)
where σ is a bandwidth parameter. To calculate the similarity w p ij in the people network G p , the demographic information, such as age, gender, address, occupation, and criminal information, such as the history of criminal records, are considered as input features. For the w c ij of the case network G c , crime reports including the location of an incident are used to measure the similarity between crime cases. In particular, to calculate the similarity between text-type crime reports, each report was converted into a term vector through term frequency-inverse document frequency (TF-IDF) [START_REF] Sanderson | prabhakar raghavan, hinrich schütze, introduction to information retrieval[END_REF] (more details are explained in Section 3). Te similarity between a person and a case is denoted as w pc ij where x i ∈ V p and x j ∈ V c . If a person is involved in a certain criminal case, then the person is connected to the case. Unlike other similarity weights in the network, the edge weight w pc ij is set to "1" (connected) if the person is involved in the case as a suspect, a victim, or a witness and "0" (disconnected) if the person is not involved in the case. Te left side of Figure 1 shows the schematic picture of the criminal network.

Crime Data.

Te crime data from January 2000 to December 2019 are collected from the Korea Information System of Criminal Justice Services (KICS) with the help of the Police Science Institute, Asan, Korea [START_REF] Kics | Korea information system of criminal justice Services (KICS)[END_REF]. Te types of crimes include a battery, assault, drug, theft, trafc violation, disorderly conduct, and fnancial crime. Te data can be categorized into people data and crime case data. Te data of people contain personal information about the suspects, victims, and witnesses. In addition, the information of past criminal history and possession of frearms appears if they exist. Te crime case data contains the date, location, type of crime, and the case summary report written by the ofcer in charge. Usually, a suspect, a victim, and a witness are related to a single case. Some cases have multiple victims or witnesses and also there are cases where a suspect from one case appears in another.

Among the crime case data variables, the summary report is unstructured data in text format. Since most of the information about the incident is included in this report, preprocessing of the text is vital. In order to process text data, nouns and verbs were extracted using the Korean morphology analyzer KoNLPy [START_REF] Park | KoNLPy: Korean natural language processing in Python[END_REF]. It is known as one of the best analyzers among open-source Korean morphology analyzers. Extracted words from the reports are converted into vectors through TF-IDF [START_REF] Sanderson | prabhakar raghavan, hinrich schütze, introduction to information retrieval[END_REF]. By using TF-IDF, the infuence of unnecessarily repeated words can be reduced to some extent, and important information can be highlighted. Figure 1 shows an example of a part of the case report, and more details are described in Table 1.

We built a criminal network with a network of 43,603 people and 20,500 cases. Both networks are sparsely connected, that is, the network densities are 2.49% and 2.90%, respectively. To overcome spare connections, the aforementioned link prediction algorithms can be used to infer missing links among suspects, victims, or witnesses to reinforce the criminal network [START_REF] Berlusconi | Link prediction in criminal networks: a tool for criminal intelligence analysis[END_REF][START_REF] Calderoni | Robust link prediction in criminal networks: a case study of the Sicilian Mafa[END_REF].

Suspect Scoring from the Criminal Network.

Te primary objective of a criminal network is to predict potential suspect candidates when a new crime case emerges, that is, identifying individuals who are likely to be involved. Suspect scores are calculated for each node in the network of people, and those with the highest scores are recommended as potential candidates. However, in the immediate aftermath of a criminal case, it is often the case where only a few individuals, such as the victim or a small number of people are directly involved, have the knowledge of the incident. Consequently, the available labeled data for training the predictive model is extremely limited, resulting in a sparse dataset. Terefore, GSSL [START_REF] Zhou | Learning with local and global consistency[END_REF] was adopted since it can learn even with just one label and further developed to ft the hierarchical structure of the criminal network.

Given a weight matrix

W � w ij 􏽮 􏽯 n i,j�1 of a plain (not layered) network G, the Laplacian is defned as L � D -W where diagonal matrix D � diag(d i ) and d i � 􏽐 n j�1 w ij .
Ten, simple GSSL works by optimizing the objective function are defned as

min f (f -y) T (f -y) + μf T Lf. (2)
Equation ( 2) minimizes the loss between predicted labels f � (f 1 , . . . , f l , f l+1 , . . . , f l+u ) T and node labels y � (y 1 , . . . , y l , 0, . . . , 0) T , while smoothing the neighbor nodes' labels to be similar, i.e., GSSL propagates the node labels to unlabeled nodes through weighted edges in the network. Te solution of (2) can be obtained as

f � (I + μL) -1 y, ( 3
)
where the parameter μ ≥ 0 controls the tradeof between loss and smoothness. Tis learning framework for a plain network is straightforward to extend to a layered network. Suppose the layered network has a total of N � n p + n c nodes and a weight matrix of W consisting of W p , W c , and W pc as shown on the right side of Figure 1, then by simply substituting the Laplacian of (2) for the plain network to that of the layered network, the solution is obtained by (3) for the layered network. If there are many layers, i.e., more than two layers are structured hierarchically, then calculating (3) is computationally highly demanding because the weight matrix is huge. In this case, approximation methods such as the Woodbury formula or the Nystrom method can be used to speed up the algorithm [START_REF] Kim | Semi-supervised learning for hierarchically structured networks[END_REF]. Te algorithm is applied to predict the suspect candidates from a small set of criminal networks [START_REF] Jhee | Fast prediction for suspect candidates from criminal networks[END_REF]. Now if case i is known, we can calculate the suspect scores. By setting the label of the i th case node as y � (0, . . . , 0, y i � 1, 0, . . . , 0) T , the label infuence propagates to both the case network and the people network, eventually computing (3). People with high f scores are regarded as the suspect candidates. Te process of suspect scoring is shown in Figure 2. First, when a new case marked in red in Figure 2(a) comes in, it is connected to the most similar existing cases. Assigning the label "1" to the new case (zeros to the remaining nodes) propagates the label to similar cases in the vicinity, and it spreads to people involved in these similar cases through the blue edge that bridges the case with people in the upper layer (Figure 2(b)). Finally, people are ranked in descending order of score (Figure 2(c)). Te questionable suspects appear as those with the highest scores (circled in red).

Fast Criminal Network-Based Suspect Prediction

Scoring suspect candidates using GSSL works well if the network is reasonably sized. However, it slows down as the number of crimes and the number of people increases, making it more likely to be lethal at a crime scene where urgent predictions are needed to solve the case. Te slowdown is mainly attributed to memory and time-consuming computations when retrieving similar cases. When a new case comes in, an exhaustive search of existing cases is conducted to fnd the most similar case. To make matters worse, every case is a high-dimensional text vector extracted from criminal case reports. For reference, the number of cases is over 20,000, the number of related people is over 100,000, and the dimensionality of a text vector is over 2,000. Ten, the search is narrowed down to only a small set of case nodes that belong to the same cluster to which the new case belongs by referencing the latent network. Te process is depicted in Figure 3(a). Terefore, the role of the latent network is important. In other words, clustering results matter. However, real data are not uniformly distributed or spherically shaped such as a normal distribution, as most well-known clustering techniques assume. Figure 3(b) illustrates clusters with diferent shapes and varying densities.

As a way to overcome this limitation, we suggest measuring the relative location of a data point by looking around its neighbors and their cluster memberships. Tat is, if the cluster memberships of its neighbors are homogeneous (belong to the same cluster), then the data point is likely to be at the core of the cluster. Conversely, if the cluster memberships of its neighbors are heterogeneous (some belong to one cluster and others belong to another cluster), then the data point is likely to be located on the boundaries of clusters. Figure 3(b) exemplifes the estimation of the relative location of data points. For x 1 , the neighbors are homogeneous in terms of cluster membership, so it is estimated to be located close to the core of cluster 1, whereas x 2 or x 3 is heterogeneous in cluster membership of its neighbors, so they are estimated to be near the cluster boundaries. Tis concept allows the cluster regions to be well-defned even when the data distribution is not spherical and the densities change.

Latent Dimension via Neighbor Mutual Information.

Te clusters are reused to reduce a higher dimension to a lower dimension. Hereafter, the resulting dimension is denoted as a latent dimension. Te size of the latent dimension is determined by the number of clusters, and the value of each dimension is determined by the degree to which the data point belongs to each cluster. Terefore, d-dimensional vectors are reduced to m-dimensional vectors that are much smaller than the original dimension (d ≫ m), and the latent values are measured by neighbor mutual information (NMI), which has been proposed here.

Te NMI of a data point (node) is the amount of mutual information between the clusters and the neighbors of the data point. Generally, for a pair of discrete random variables X and Y, mutual information is defned as 􏽐 x∈X 􏽐 y∈Y p(x, y)logp(x, y)/p(x)p(y), where p(x, y) is the joint probability [START_REF] Kraskov | Estimating mutual information[END_REF][START_REF] Bachman | Learning representations by maximizing mutual information across views[END_REF]. If p(x, y) � p(x)p(y), then the logarithm becomes zero and thus mutual information becomes zero, which means X and Y are independent. Analogously, the NMI of node x i is defned between the cluster C m ∈ C 1 , . . . , C M 􏼈 􏼉 and the set of its neighbors Ne(x i ) as

NMI x i , C m 􏼁 � Ne x i 􏼁 ∩ C m 􏼌 􏼌
| � N), then NMI(x i , C m ) � N • |Ne(x i )|/|Ne(x i )| • |C m | � 1,
whereas when none of the neighbor nodes belong to the cluster, then, NMI(x i , C m ) � 0. Finally, z ij represents the normalized version of NMI(x i , C m ), which satisfes the nonnegativity and sum to one constraint. Tere are tradeofs between the number of clusters and the latent dimension. If the number of clusters increases, the number of centroids in the latent network increases, and it leads to a reduction in the scope of local search because the data are split into more clusters and thus fewer data in each cluster. However, the latent dimension increases, that is, increasing the number of clusters afects both increasing and decreasing computation time and vice versa. Terefore, it is necessary to adjust the number of clusters according to the situation at hand, so as to fgure out whether search coverage or dimensionality reduction is more important.

z im � NMI x i , C m 􏼁 􏽐 m NMI x i , C m 􏼁 . ( 5
)

Scoring via the Latent Network. Te latent network makes the computation of node scoring for the case network much lighter. Given a weight matrix 􏽢

W � 􏽢 w ij 􏽮 􏽯 M i,j�1 between M centroids, a score vector 􏽢 f � ( 􏽢 f 1 , . . . , 􏽢 f M )
T can be calculated similarly to [START_REF] Yan | Research on prediction of multi-class theft crimes by an optimized decomposition and fusion method based on XGBoost[END_REF]. It is now straightforward to derive a score vector f i � z T i 􏽢 f, which is the sum of the centroid scores 􏽢 f weighted by the latent vector z i . Terefore, the scoring problems (2) and (3) turn into fnding the optimal centroid scores on the latent network. Te objective function is represented as follows:

min 􏽢 f (Z 􏽢 f -y) T (Z 􏽢 f -y) + μ 􏽢 f T 􏽢 L 􏽢 f, ( 6
)
where 􏽢 L � Z T LZ ∈ R M×M , i.e., the graph Laplacian of the latent network. Te only diference from ( 2) is that the loss is computed using the weighted sum of centroid scores and the smoothness is optimized on the latent network. It is thus trivial to derive the solution to [START_REF] Vijayarani | A comprehensive analysis of crime analysis using data mining techniques[END_REF].

􏽢 f � Z T Z + μ 􏽢 L 􏼐 􏼑 -1 Z T y. ( 7
)
Finally, the score vector f is predicted using the centroid scores 􏽢 f and the latent representation, Z, that is,

f � Z 􏽢 f 1 T Z 􏽢 f , ( 8
)
where 1 T Z 􏽢 f is the normalization. From now on, the proposed method is named as MISSL, that is, the abbreviation of neighborhood information-based semisupervised learning. Constructing the latent network takes an extra O(NM 2 ), but since M is small, it is not a huge burden. However, once the latent network is constructed, it provides signifcant advantages for computing the inverse matrix of the solution. For the latent network in [START_REF] Sparrow | Te application of network analysis to criminal intelligence: an assessment of the prospects[END_REF], it has O(M 3 ) complexity, whereas the original network in (3) has O(N 3 ) complexity (N ≫ M). Also, MISSL was originally designed to work on undirected networks. However, it can be extended to directed networks by converting asymmetric matrices to the symmetric graph Laplacian [START_REF] Zhou | Semi-supervised learning on directed graphs[END_REF].

Application to Criminal Network.

Applying MISSL to criminal networks is simple. Conceptually, a latent network of centroids lies between the people network and the case network, connecting the two networks. Edges from a centroid node to case nodes are connected via the latent vector, as explained in the previous section. When a new case comes in, MISSL fnds the centroid of the nearest cluster C m , and converts node x new into a new latent vector using only the representation matrix for that cluster. More precisely, we frst calculate the similarity between centroids and the node x 1 is estimated to be near the core of cluster 1 because all three of its neighbors belong to cluster 1 (homogeneous in terms of cluster membership), whereas x 2 or x 3 is estimated to be near the boundaries of clusters because the neighbors difer in cluster membership (heterogeneous). (c) Te cluster memberships of x 1 , x 2, and x 3 , NMI calculation result, and their latent vector Z.

x 1 x 2 x 3 C 1 C 2 C 3 |Ne (x i )| 3 2 1 0 1 1 0 0 1 x 1 x 2 x 3 C 1 C 2 C 3 (c)
x new to fnd the nearest cluster. Within the cluster C m , mixing weight ω is optimized to transform x new to z new and the latent vector Z is updated. Ten, we applied MISSL for the new case using the updated latent vector. Conversion from x new to z new is performed by fnding a mixing weight ω ∈ R M×d which is obtained by solving the following problem:

min ω Z m ω -X m 􏼁 T Z m ω -X m 􏼁, (9)
where X m ∈ R N m ×d and Z m ∈ R N m ×M are the matrices of the data points belonging to cluster C m . From this, the optimal weight for conversion is calculated as ω � (Z T m Z m ) -1 Z T m X m . So, the latent vector for the new case is

z new � ωx new � Z T m Z m 􏼐 􏼑 -1 Z T m X m 􏼔 􏼕x new . ( 10
)
Te new latent vector is connected to the case network by fnding the closest cases. MISSL then scores suspect candidates in the people network according to [START_REF] Sparrow | Te application of network analysis to criminal intelligence: an assessment of the prospects[END_REF]. By using MISSL, a much faster prediction is possible than the naïve SSL in Section 2. Tis allows ofcial crime investigators to quickly weed out suspect candidates at a crime scene.

Experiments on Benchmark Data

Experimental results in the following section show that the proposed algorithm has a fast inference time and competitive performance compared to the existing approaches.

4.1. Data. We evaluated the proposed method, MISSL, on benchmark datasets. Te datasets are g241c, g241n, MNIST, and CIFAR-10. g241c and g241n datasets were artifcially generated to hold the cluster assumption. Te data points of g241c were drawn from each of the two unit-variance isotropic Gaussians. Te label of a data point represents the Gaussian it was drawn from. Te data points of g241n were drawn from each of the two unit-variance isotropic Gaussians which have a potentially misleading cluster structure and no manifold structure. Te centers for the positive class have a distance of six in a random direction and the centers for the negative class were fxed by moving from the former centers to a distance of 2.5 in a random direction. All dimensions were standardized to zero mean and unit variance. Te number of dimensions and data points are 241 and 1,500, respectively [START_REF] Chapelle | Semi-supervised learning[END_REF]. MNIST is a handwritten digit dataset of 28 by 28 grayscale normalized and centered images. Te labels of the dataset contain the number from zero to nine and the number of data points is 70,000 [START_REF]MNIST handwritten digit database, yann LeCun, corinna cortes and chris burges[END_REF]. Te CIFAR-10 dataset (Canadian Institute for Advanced Research) is a subset of the 80 million tiny images. Te CIFAR-10 dataset contains 60,000 32 × 32 color images in 10 diferent classes. Te 10 diferent classes represent airplanes, cars, birds, cats, deer, dogs, frogs, horses, ships, and trucks. Tere are 6,000 images of each class [START_REF] Krizhevsky | Learning multiple layers of features from tiny images[END_REF]. Table 2 summarizes the datasets.

Experimental Setup.

We compared the performance and scalability of MISSL with another scalable method called anchor graph regularization (AGR) [START_REF] Liu | Large graph construction for scalable semi-supervised learning[END_REF]. Here, an anchor of AGR that employs RBF kernel or local linear embedding corresponds to a centroid of our method [START_REF] Roweis | Nonlinear dimensionality reduction by locally linear embedding[END_REF]. Te two network-based methods use anchors or centroids to reduce the heavy computation of GSSL. Te performance of Naïve GSSL was used as a reference. Tis indicates that both methods, MISSL and AGR, cannot overwhelm the baseline performance. Te dataset was divided into the 20-cross validation (20-cv) format. However, unlike the supervised learning setting, SSL assumes a few labeled data, so a single set is the training set, and the remaining 19 sets are confgured as the validation sets; that is, the labeled data is 5% (positive: 2.5% and negative: 2.5%) of the entire dataset. Te experiment was repeated 10 times to optimize the hyperparameters. 3 shows the performance results on benchmark datasets measured with the area under the receiver operating characteristic curve (AUC), along with the hyperparameters of AGR, m (the number of centroids) and r (the number of nearest centroids), and hyperparameters of MISSL, m (the number of clusters) and k (the number of nearest neighbors). Te average performances over 20-cv of AGR and MISSL are 97.65% and 98.98% of the reference performance, respectively. MISSL showed better performance than AGR in all datasets. It shows that the latent dimension using mutual information of MISSL better represents data than RBF kernels or LLE of AGR.

Results. Table

Te experimental results on computation time are shown in Table 4. Te computation time was measured separately in three parts: the network construction (or clustering), the latent representation, and the inference (prediction). For GSSL, no conversion time for new representation is required, and clustering time is only required for AGR and MISSL. However, GSSL takes a lot of time for both network construction and inference, whereas AGR or MISSL signifcantly reduces computation time by using clustering. Compared to AGR, MISSL is inherently superior to AGR in terms of speed as it does not require an optimization process. A comparison of the overall time empirically proves that MISSL is more efcient. For the small datasets, such as g241c and g241n, MISSL and AGR are 100 times and 90 times faster than the reference time, respectively, and for the large datasets, MNIST and CIFAR-10 are 1,000 times and 900 times faster.

Applications for Suspect Candidate Scoring Using the Criminal Network

We applied MISSL to real crime data provided by the Korean National Police. Te experimental setup and results are as follows.

Criminal

Results of Suspect Candidate Scoring.

In Figure 5, the AUCs for 100 cases are shown. Te overall average AUC for 20,500 cases is 0.82 ± 0.01. More specifcally, Figure 6 presents a typical scoring curve. Suspect candidates for case #8 are sorted along the curve. Te blue nodes represent predicted suspect candidates and the red represents the actual suspect. As shown in fgure, the real suspect is ranked relatively high, within the top 1%. Case #12 was similar to case #8. Both cases are similar in that the suspect steals money or jewelry. Te summary reports of cases #8 and #12 are described on the right side with the personal information of the suspect. Indeed, the suspect of case #12 was highly scored (red circle) and was the actual ofender of case #8. Figure 7 shows the computation time and performance of MISSL while varying the number of clusters or centroids, respectively. GSSL are indicated by red and blue dotted lines in the fgure. By comparing the computation time, it was found that the MISSL was 0.06 seconds when the number of clusters was 100, which is about 1.76e + 04 times faster than that of GSSL (1,058.65 seconds). Te computation time increases as the number of clusters increases but is still very trivial compared to GSSL. Te performance of MISSL increases as the number of clusters increases. Te highest AUC of 0.8 was reached when the number of clusters was 100. From the perspectives of both computational time and performance, MISSL provides a high accuracy within a reasonable amount of time. Tis is critical in real-world applications. In an urgent case, the police do not have to wait all day to get a list of suspect candidates.

Conclusion and Discussion

In this study, we proposed a framework for predicting suspect candidates based on the criminal network. Te algorithm we employed is graph-based SSL, which may be inappropriate when networks are large and complicatedly structured. So, to put the GSSL to practical use in the criminal network, we developed an algorithm based on latent representation and mutual information. Te proposed method, MISSL, shows almost similar performance to the graph-based SSL but has a much faster inference time. As an application, a criminal network is constructed from realworld crime data, and suspect candidate scoring is performed by MISSL. Te predicted results show the validity and efcacy of MISSL. Te framework of suspect candidate scoring introduces a novel way of analyzing crime data, and the results shed a new light on network-based machine learning approaches for social network analysis. Future eforts may identify a more efcient mechanism to optimize the hyperparameter of MISSL and the number of clusters (i.e., the number of latent dimensions) that afect performance. Empirically, a higher AUC was obtained from the increased number of clusters. From another perspective, clustering can be expensive for large-scale datasets. Terefore, further research on faster clustering methods or modeling the algorithms robust to clustering will enrich the results of our study.

p 1 p 36 Figure 1 :

 1361 Figure 1: Te schematic picture of the criminal network. Te upper layer is the people network and the lower layer is the case network. Te corresponding weight matrices are represented as W p and W c in the fgure on the right, respectively. Two networks are connected by person-to-case edges denoted by W pc in the weight matrix. A person node has demographic information of an individual while a case node has information on crime reports (texts are translated from Korean to English. Names or titles in texts are de-identifed. All the meaningless symbols were removed).

Figures 3 ("Figure 2 :

 32 Figures 3(b) and 3(c) show several examples of NMIs that transform x i to the relative locations between clusters. For node x 1 , NMI(x 1 , C 1 ) � (30 × 3)/(3 × 12) � 2.5 by (4) since C 1 has 12 nodes and all three neighbors are a part of the cluster, whereas NMI(x 1 , C 2 ) and NMI(x 1 , C 3 ) are zeros. Tus, z 1 � [1, 0, 0]. Note that, node x 1 is in the core of C 1 far from C 2 and C 3 . Conversely, node x 3 is located on the boundary of C 2 and has neighbors belonging to C 1 , C 2 , and

Figure 3 :

 3 Figure 3: Te latent network and the latent vectors. (a) A criminal network with the latent network composed of centroids of the case clusters. Colors represent cluster membership and squares represent cluster centroids. (b) Te relative location estimation in circumstances of clusters of diferent shapes and densities.x 1 is estimated to be near the core of cluster 1 because all three of its neighbors belong to cluster 1 (homogeneous in terms of cluster membership), whereas x 2 or x 3 is estimated to be near the boundaries of clusters because the neighbors difer in cluster membership (heterogeneous). (c) Te cluster memberships of x 1 , x 2, and x 3 , NMI calculation result, and their latent vector Z.

Figure 4 :

 4 Figure 4: Te crime case network of 20,500 cases. Te node color represents the type of crime. Less than 1% of types are grouped by "etc." Te size of a node indicates the node degree. Te subnetwork of "assault" is depicted on the right insert. Te number in the case node represents the case report number. Te cases extracted show a similar criminal pattern in which the suspects either beat or detain the victims.

Figure 7 :

 7 Figure 7: Te computation time (line chart) and AUC performance (bar chart) of MISSL while varying the number of clusters. Te performance of GSSL is shown by red dotted lines.

Table 1 :

 1 Crime data description.

	Variables

Type Battery, assault, drug, theft, trafc violation, disorderly conduct, and fnancial crime similarities to the centroids and selects the closest cluster.

Table 2 :

 2 Benchmark datasets.

	Dataset	Classes	Dimension	Nodes	Remark
	g241c	2	241	1,500	Artifcial
	g241n	2	241	1,500	Artifcial
	MNIST	10	784	70,000	Large-scale
	CIFAR-10	10	3,072	60,000	Large-scale

Table 3 :

 3 AUC comparison on benchmark datasets.

	Method	g241c	g241n	MNIST	CIFAR-10
	AGR	0.6438 ± 0.0288	0.6448 ± 0.0468	0.9876 ± 0.0002	0.7305 ± 0.0014
	(m, r)	(150, 5)	(150, 5)	(300, 10)	(300, 10)
	MISSL	0.7652 ± 0.0266	0.6744 ± 0.0257	0.9885 ± 0.0006	0.7324 ± 0.0022
	(m, k)	(150, 10)	(150, 10)	(300, 20)	(300, 20)
	GSSL	0.7808 ± 0.01890	0.7504 ± 0.0682	0.9896 ± 0.0003	0.7340 ± 0.0025
	Te bold values refer to the values of the proposed method (MISSL).

Table 4 :

 4 Computation time (sec.) on benchmark datasets.

	Data	Method	Network construction (or clustering)	Latent representation Z	Inference	Total
		AGR	0.20	1.19	0.01	1.40
	g241c	MISSL	0.20	0.08	0.01	0.29
		GSSL	0.25	-	0.06	0.31
		AGR	0.21	1.19	0.01	1.41
	g241n	MISSL	0.21	0.08	0.01	0.30
		GSSL	0.25	-	0.06	0.31
		AGR	25.00	67.29	0.08	92.37
	MNIST	MISSL	25.00	0.09	0.07	25.16
		GSSL	482.48	-	175.71	658.19
		AGR	85.12	77.70	0.02	162.84
	CIFAR-10	MISSL	85.12	0.06	0.01	85.19
		GSSL	1,221.03	-	111.54	1,332.57

Te bold values refer to the values of the proposed method (MISSL).

  Network and Experimental Setting. Figure4visualizes clusters in the case network using OpenOrd[START_REF] Martin | OpenOrd: an open-source toolbox for large graph layout[END_REF]. Te node color represents the type of crime. Te same type of crime is grouped closely in the network and batteries and assaults are grouped together since they are both violent crimes. For validation on suspect scoring, a case node is randomly selected as a test node, and the edges connected to the node are removed from the network (including the edge to the actual suspect). Te top 20 people with the highest scores were reported as suspect candidates. Te performance is measured by calculating the AUC of the predicted suspects for each crime case. Te leave-one-out (LOO) method is applied to 20,500 cases. AUC of suspect candidate scoring for each crime case. 100 cases are randomly sampled from 20,500 cases. Te worst case shows around 0.65 AUC and the best case shows around 0.93 AUC. Scores of suspect candidates for case #8. Blue circles indicate suspect candidates, and a red circle indicates actual suspects. Te suspect is in the top 1% of suspect candidates. Case #12 was similar to case #8.

	1
	0.9
	0.8
	AUC
	0.7
	0.6
	0.1
	0
	244	653	707	733	947	1106	1556	1992	2000	2440	2604	2664	2842	2909	3061	3232	3334	3510	3831	4031	4589	4993	5148	5213	5230	5279	5659	5677	5710	5860	6501	6912	6978	7175	7210	7800	7822	8041	8647	8995	9135	9623	9703	9951	10041	10217	10373	10882	11212	11218	11270	11641	11661	11998	11998	12629	12964	13250	13430	13438	13443	13915	13934	14244	14332	14474	14542	15235	15401	15452	15472	15523	15534	15694	15973	16241	16302	16406	16693	16702	16881	17032	17235	17408	18264	18569	18725	18773	18803	19050	19147	19148	19480	19622	19629	19666	19670	19675	19781	19898
																																																Crime case
	Figure 5: 1.0																																									Suspect
					0.6
	Score				0.5 0.4
																																																				1%
					0.3
					0.2
					0.1
					0.0
					Suspect
					Suspect candidate
	Figure 6:

International Journal of Intelligent Systems

International Journal of Intelligent Systems

Acknowledgments

Tis research was supported by the Institute for Information Communications Technology Promotion (IITP) grant funded by the Korea Government (MSIT) ((Grant no. 2022-0-00653) Voice Phishing Information Collection and Processing and Development of a Big Data-Based Investigation Support System), BK21 FOUR program of the National

Data Availability

Te data that support the fndings of this study are available from the corresponding author upon reasonable request.

Conflicts of Interest

Te authors declare that they have no conficts of interest.

Research Foundation of Korea funded by the Ministry of Education (Grant no. NRF5199991014091), and the Ajou University research fund.


Huu-Thinh Do
email: huu-thinh.do@lcis.grenoble-inp.fr

Franco Blanchini
email: franco.blanchini@uniud.it

Stefano Miani
email: miani.stefano@uniud.it

Ionela Prodan
email: ionela.prodan@lcis.grenoble-inp.fr

LP-generated Control Lyapunov Functions with application to multicopter control

Keywords: Constrained control, Lyapunov function, open loop trajectories, linear programming, homogeneous polynomials

In this work, we study a technique of exploiting open-loop generated trajectories for a constrained control problem, using them to shape suitable non-quadratic Control Lyapunov Functions. These trajectories, generated off-line, allow detecting a suitable domain of attraction in which a candidate Lyapunov function has a negative derivative. Given a suitably constructed basis function, our working machinery is based on linear programming, hence, the technique can be applied to problems of non-trivial size in terms of the number of basis functions and points in the state space. For linear systems, we seek convex Lyapunov functions, which are homogeneous polynomials. Simulation and experimental results for drone control are given.

I. INTRODUCTION

The Control Lyapunov Function (CLF) notion has proven fundamental in control design, especially when uncertainty and constraints come into play. Determining a CLF is a pivotal step, after which a controller can be immediately deduced using Sontag's universal formula [START_REF] Sontag | A 'universal'construction of Artstein's theorem on nonlinear stabilization[END_REF]. In view of applications with input constraints, several explicit control laws are proposed, depending on the constraint's characterization (e.g. Minkowski balls [START_REF] Malisoff | Universal formulas for feedback stabilization with respect to minkowski balls[END_REF], their direct products [START_REF] Kidane | Controller design for a nonlinear system with inputs restricted to a direct product of minkowski balls[END_REF] or recently, polytopes [START_REF] Yamashita | Asymptotic stabilization of nonlinear systems with convex-polytope input constraints by continuous input[END_REF], [START_REF] Leyva | Stabilization of affine systems with polytopic control value sets[END_REF]). Furthermore, via a Quadratic or Linear Program, implicit continuous control can be introduced with different objectives (min-norm [START_REF] Freeman | Inverse optimality in robust stabilization[END_REF] or max-rate [START_REF] Suarez | Global clf stabilization for systems with compact convex control value sets[END_REF]

control).

However, proposing a CLF a priori to the controller synthesis is not a trivial task, especially when its stabilizing condition (negative time derivative) has to be proven feasible under operating constraints. This task is usually handled on a case-by-case basis, with specific numerical tools and some restrictions on the number of states and inputs. For example, with sum-of-square methods, CLF for polynomial systems can be achieved with some relaxation via semidefinite programming (SDP) [START_REF] Prajna | Nonlinear control synthesis by sum of squares optimization: A lyapunov-based approach[END_REF], [START_REF] Chesi | LMI techniques for optimization over polynomials in control: a survey[END_REF]. In [START_REF] Hu | Composite quadratic lyapunov functions for constrained control systems[END_REF], a composite quadratic function was developed as the CLF for linear systems subject to saturation and state constraints, resulting in the domain of attraction described by a convex hull of ellipsoids. The study for linear piecewise affine systems is conducted in [START_REF] Lazar | On infinity norms as lyapunov functions for piecewise affine systems[END_REF]. More generally for nonlinear systems, yet at the price of high computational burden with respect to the dimensionality problem, procedures to numerically compute a continuous piecewise affine CLF were developed in [START_REF] Lavaei | Systematic, lyapunov-based, safe and stabilizing controller synthesis for constrained nonlinear systems[END_REF], [START_REF] Steentjes | Construction of continuous and piecewise affine feedback stabilizers for nonlinear systems[END_REF] via state space triangulation. In [START_REF] Min | Data-driven control with inherent lyapunov stability[END_REF], under a known assumed approximation error bound, a neural network structure was proposed to learn the dynamics, the Lyapunov function, and the controller simultaneously. Online synthesis of state-dependent CLF as well as the corresponding control is proposed for discrete systems [START_REF] Lazar | Synthesis of trajectory-dependent control Lyapunov functions by a single linear program[END_REF], [START_REF] Lazar | Flexible control lyapunov functions[END_REF]. For certain low-dimensional systems, the Hamilton-Jacobi-Bellman partial differential equation can be solved to incorporate the notion of optimality [START_REF] Beard | Galerkin approximations of the generalized Hamilton-Jacobi-Bellman equation[END_REF], [START_REF] Aguilar | Numerical solutions to the Bellman equation of optimal control[END_REF]. For a more detailed survey of computational CLF synthesis methods, we refer to the review [START_REF] Giesl | Review on computational methods for lyapunov functions[END_REF].

Meanwhile, one effective strategy to sidestep the a priori design of either the controller in closed form or a CLF is the Model Predictive Control (MPC) approach. In the standard setting, for example, a proper choice of terminal ingredients (cost and constraints) will result in a stabilizing implicit controller, and a cost function that can be interpreted as a piecewise CLF [START_REF] Bayer | A tube-based approach to nonlinear explicit mpc[END_REF]. But then, the real-time implementation will become more time-consuming due to the solving of an online optimization problem. Furthermore, for controllers with limited hardware power, this difficulty may lead to computational unfeasibility or instability due to the so-called early termination [START_REF] Hosseinzadeh | Robust to early termination model predictive control[END_REF] (i.e., the search for the optimal control is interrupted due to the limited execution time).

With those pros and cons of both CLF-based control and MPC strategy, our objective is to exploit the systematically designed MPC law to seek a corresponding CLF during the offline synthesis. More specifically, in this work, we propose a synthesis procedure for a CLF by parameterizing it as linear combinations of certain candidate CLFs. The idea includes collecting admissible state-input pairs generated by a stabilizing MPC software, then a parameterization is constructed by a linear program. The advantages of the technique are the simplicity of the offline synthesis and the reduction in computational overhead during the online implementation, allowing applications to high-dimensional systems. Then, to validate its applicability and emphasize the computational advantages, the method is experimentally examined over a challenging and timely application, the quadcopter position stabilization. It is well known that it requires small sampling time, online feasibility, and constraint satisfaction. The computational requirements become even more stringent when considering a team of drones. Summarily, in this work, we:

• propose an offline procedure to construct a CLF based on linear programming, exploiting the admissible control generated by an MPC law; • validate the proposed process for a drone stabilizing problem through experimental tests. The video for the experiment can be found at https://youtu.be/ XPyZAJoN2Ps.

The remainder of the paper is structured as follows. Section II presents the preliminaries and the main principles for the CLF construction. Therein, via a linear program (LP), a CLF is constructed with stabilizing conditions imposed at points in the state-input space. Section III presents the formulation for the acquisition of admissible points shaping the desired CLF from an MPC-based admissible point generator. Section IV particularizes the proposed technique for a quadcopter position control problem. The resulted CLF and the associated control will be validated in Section V with various experimental tests. Finally, Section VI concludes and discusses future directions.

Notation: Bold lower-case letters denote vectors. With

x = [x 1 , x 2 , ..., x n ] ⊤ ∈ R n , ∇V (x) = [ ∂V (x)
∂x1 , ∂V (x) ∂x2 , ..., ∂V (x) ∂xn ] denotes the gradient of the value function V (x). The standard norms for a vector x are denoted as ∥x∥ 2 ≜

√

x ⊤ x, ∥x∥ P ≜ √

x ⊤ P x. ∥ • ∥ ∞ is the infinity norm. 0 denotes the origin of the vector space with the appropriate dimension. The letter k denotes the discrete step, while t represents the time variable. N (n, a, b) represents a finite set of n evenly spaced real numbers over the interval [a, b]. co{•} denotes the convex hull. diag(•) returns a matrix with its entries diagonally arranged.

II. MAIN IDEA AND MACHINERY

Consider the following nonlinear control affine system:

ẋ = f (x) + g(x)u, (1)
with a convex and compact constraint set:

u ∈ U, (2)
where x ∈ R n denotes the state vector, and u ∈ R m is the input vector constrained in U.

Assumption 1: Assume the pair x = 0 and u = 0 is a steady state, i.e., U includes 0 in its interior and f (0) = 0. □ Assume that we have an open-loop constrained trajectory software tool, capable of driving the system to 0. Then we could just use the software in an MPC fashion for control. Still, for systems with a small sampling time or limited hardware resources, such MPC implementation could be challenging.

The underlying idea is that of adopting this open-loop trajectory generator to simply determine candidate state controlpair optimal trajectories. These state control pairs are then used to generate suitable constraints for finding a CLF, within a linearly parameterized family, through linear programming. These generated trajectories will eventually be disregarded, and only the CLF will be retained for actual control. To this aim, some standard definitions are recalled in the following.

Definition 1 (Control Lyapunov Function): A positive definite smooth Lyapunov function, V (x), is a CLF for (1) if there exists a function u = Ξ(x) with values in U and dV (x)/dt = ∇V (x)[f (x) + g(x)Ξ(x)] ≤ -βV (x), (3) for some β > 0.

□ Definition 2: The state-control pair (x, u) for (1), u ∈ U is admissible for V (x) if it satisfies [START_REF] Kidane | Controller design for a nonlinear system with inputs restricted to a direct product of minkowski balls[END_REF]. □ Definition 3 (Domain of Attraction): D ⊆ R n is a domain of attraction (DOA) of system (1) if there exist u(t) = Ξ(x(t)), such that any trajectory x(t) with x(0) ∈ D, satisfies:

x(t) ∈ D, u(t) ∈ U, and lim

t→∞ x(t) = 0. (4)

□

Next, in the offline control synthesis, adopting the available software, we generate several admissible pairs and then, formulate the CLF based on these points. This is a simple idea to generate, under appropriate conditions, points that are inside the domain of attraction D, eliminating all points outside, for which, no constrained control can bring the state to 0.

Let us assume that the open-loop trajectory generator ensures the existence of constrained trajectories x(t), u(t), from which we can sample M pairs (x j , u j ), u j ∈ U, j ∈ {1, ..., M }. Consider the following candidate CLF:

V (x) = N i=0 α i V i (x), α i ≥ 0, (5)
where V i (x) are smooth positive definite Lyapunov functions. Next, consider the following linear programming problem with N variables, and M inequality constraints formulated from M state-input pairs collected previously:

min αi N i=1 ℓ i α i (6a) N i=1 α i [∇V i (x j ) [f (x j ) + g(x j )u j ] + βV i (x j )] ≤ 0, (6b) N i=1 α i = 1, α i ≥ 0, j ∈ {1, ..., M }, (6c)
where ℓ i are weights. These weights give freedom of choice in the computation. As an example, if V 1 (x) is a "privileged" function, for instance a quadratic function associated with a local control optimizing cost, then ℓ 1 can be adjusted to influence the function dominance within the parameterization [START_REF] Leyva | Stabilization of affine systems with polytopic control value sets[END_REF]. Then, as a result of (6), we achieve a parameterized CLF established as a linear combination of V i (x) with the coefficients α i . An illustration on how the Lyapunov function V (x) is shaped by the stabilizing constraint (6b) at each sampled data pairs is given in Fig. 1. Remark 1: With the intuition of parameterizing the function V (x) via the set of V i (x), it is essential to provide a suitable choice for such candidates. This, indeed, can be done via an iterative process where, after one solution of (6), the function V i (x) with α i ≈ 0 can be discarded and a new set of functions can be adapted. Note that, this can be accomplished thanks to the simplicity of the LP [START_REF] Freeman | Inverse optimality in robust stabilization[END_REF]. A discussion on the choice of V i (x) will be given later in the illustrative example.

□ Forming the CLF in this manner has some pros and cons:

• by employing an LP, it is possible to use a considerably large number of points and basis functions V i (x), ensuring the existence of α i for some β in [START_REF] Freeman | Inverse optimality in robust stabilization[END_REF], even in high-dimensional systems; • although the number of generating data pairs (x j , u j ) can be huge, we are still working on a finite set.

To address a preliminary answer to the second point, let us consider the following propositions. For the set U, denote by Ψ U (w) its support functional:

Ψ U (w) ≜ max u∈U w ⊤ u. ( 7
) X V (x) ∇V (x j ) ⊤ f (xj) + g(xj)uj (xj, uj) CLF constraint (6b) ∇V (xj)(f (xj) + g(xj)uj) ≤ -βV (xj) 1 Collect admissible pairs •(xj, uj)
Impose the constraints & solve the program (6) for

V (x) = N i=0 αiVi(x) 2 3 - 4
Verify the condition ( 8)

Fig. 1. Geometric interpretation of how the stabilizing condition (6b) is imposed from the sampled pairs (red points) (x j , u j ) to the surface of V (x).

Note that Ψ U (0) = 0 and that, in view of Assumption 1, Ψ U (w) is positive definite and positively homogeneous of degree one:

Ψ U (λw) = λΨ U (w), ∀ λ > 0.
Proposition 1: The function V (x) is a CLF on a given domain X containing 0 in its interior iff it satisfies the following inequality on X :

∇V (x)f (x) -Ψ U (g(x) ⊤ ∇V (x) ⊤ ) + βV (x) ≤ 0, (8)
for some β > 0 on X . □ Proof: By Definition 1, V (x) is a CLF iff ∀x ∈ X , ∃u ∈ U such that:

∇V (x)[f (x) + g(x)u] + βV (x) ≤ 0, (9)
or, equivalently:

min u∈U ∇V (x)[f (x) + g(x)u] + βV (x) ≤ 0. (10)
Condition [START_REF] Hu | Composite quadratic lyapunov functions for constrained control systems[END_REF] then can be rewritten as:

max u∈U ∇V (x)f (x) -∇V (x)g(x)u + βV (x) ≤ 0. (11)
Then, by the support function definition in [START_REF] Suarez | Global clf stabilization for systems with compact convex control value sets[END_REF], one gets to Proposition 1. Remark 2: It is worth noticing that (1) does not refer to a specific control action u = Ξ(x) but depends only on the shape of U. For instance, if, as often assumed in the literature, U is the unit ball ∥u∥ 2 ≤ 1, then equation ( 8) becomes:

∇V (x)f (x) -∥g(x) ⊤ ∇V (x) ⊤ ∥ 2 + βV (x) ≤ 0 (12)
□ For a domain of interest X , the proposition is useful since it allows testing the generated function V (x) in ( 6), by dense gridding and verifying the condition (8) in X . To recap, we propose the following procedure to determine V (x).

Procedure 1 Offline CLF generation 1: Collect M admissible pairs (x j , u j ) from a stabilizing controller for system (1). This can be sampled pairs from the system's trajectories with several random initial conditions (see also Fig. 1); 2: Fix a basis function V i (x) as in ( 5), i = 1, ..., N of positive definite functions, and β > 0 as in (3); 3: Formulate and solve the linear program (6); 4: If the problem ( 6) is feasible, then grid the state space within a region of interest and check ( 8) for all points in the grid.

Aligning with the common CLF setup, the aforementioned considerations are formalized in the following proposition.

Proposition 2: Suppose that Procedure 1 is successful and that there exists κ > 0 for which inequality ( 8) is verified in the sublevel set D κ = {x : V (x) ≤ κ}, then V (x) is a CLF for the problem (1) with the DOA D κ . □ Besides intensely gridding the state space inside D κ to verify [START_REF] Prajna | Nonlinear control synthesis by sum of squares optimization: A lyapunov-based approach[END_REF], one can also achieve an inner estimation as follows.

Proposition 3: Consider an ellipsoid inscribed in the input constraint set U in ( 2), which has the form:

E = {u = Ew, ∥w∥ 2 ≤ 1} and E ⊂ U. (13)
Let Φ(κ e ) be defined as:

Φ(κ e ) ≜ max V (x)≤κe ∇V (x)f (x) -∥Eg(x) ⊤ ∇V (x) ⊤ ∥ 2 + βV (x) (14)
then, for κ e > 0, if Φ(κ e ) ≤ 0, D κe is a DOA. □ Proof: From the convexity and compactness of E and U, one can state [START_REF] Vitale | Lp metrics for compact, convex sets[END_REF]:

E ⊂ U ⇒ Ψ E (v) ≤ Ψ U (v), ∀v ∈ R m . ( 15
)
Thus,

∇V (x)f (x) -Ψ U (g(x) ⊤ ∇V (x) ⊤ ) + βV (x) ≤ ∇V (x)f (x) -Ψ E (g(x) ⊤ ∇V (x) ⊤ ) + βV (x) = ∇V (x)f (x) -∥Eg(x) ⊤ ∇V (x) ⊤ ∥ 2 + βV (x). (16)
Consequently, if for κ e > 0, Φ(κ e ) ≤ 0, then, by Proposition 1, D κe is a DOA. With this method, we iteratively seek also for the largest κ e that satisfies Φ(κ e ) ≤ 0. The trade-off for this is the conservative employment of an ellipsoidal subset E instead of the original constraint set U.

To complete the synthesis, the remaining discussion resides on how to formulate the generator of admissible pairs (x j , u j ) for [START_REF] Freeman | Inverse optimality in robust stabilization[END_REF]. Moreover, with the established CLF, the choice of the associated control law is an important question as well. These two issues will be addressed in the next section.

III. ADMISSIBLE PAIR GENERATOR AND CONTROL

IMPLEMENTATION

In this part, the collection of admissible pairs at step 1 in Procedure 1 is first discussed. Next, the online control policy employed with the generated CLF will be introduced.

A. Admissible pair generation

As mentioned before, we aim to employ an MPC solver to achieve stabilizing trajectories from which admissible pairs are sampled. However, the MPC strategy commonly needs to be implemented in a discretized manner. With such a discretization, the generated pairs do not necessarily imply the stability property or the accurate evolution of the continuous system [START_REF] Sontag | A 'universal'construction of Artstein's theorem on nonlinear stabilization[END_REF]. For this reason, we employ the Euler Auxiliary System (EAS) [START_REF] Blanchini | Suboptimal receding horizon control for continuous-time systems[END_REF]:

x(k + 1) = x(k) + τ [f (x(k)) + g(x)u(k)] , (17)
where τ > 0 is the discretization constant parameter. Then, with this setup and a CLF generated with [START_REF] Beard | Galerkin approximations of the generalized Hamilton-Jacobi-Bellman equation[END_REF], we have the following result. Theorem 1: Any convex function V (x) which is a CLF for the EAS in ( 17) is also a CLF for the associated continuous time system in [START_REF] Sontag | A 'universal'construction of Artstein's theorem on nonlinear stabilization[END_REF] .

□ Proof: The convexity of V (x) implies that the difference quotient function is non-decreasing with regard to h ∈ R + and the Dini derivative yields [START_REF] Rockafellar | Convex analysis[END_REF]:

d dt V (x) = lim h→0 + sup V (x + h[f (x) + g(x)u]) -V (x) h ≤ V (x + τ [f (x) + g(x)u]) -V (x) τ , for τ ≥ h. (18)
Meanwhile, being a CLF of the EAS implies, ∃β d > 0:

V (x + τ [f (x) + g(x)u]) -V (x) ≤ -β d V (x). (19)
Thus, from [START_REF] Aguilar | Numerical solutions to the Bellman equation of optimal control[END_REF], [START_REF] Giesl | Review on computational methods for lyapunov functions[END_REF] and Definition 1, one can state that :

d dt V (x) ≤ - β d τ V (x), (20)
which completes the proof. Therefore, for the offline admissible pair generator in Procedure 1, step 1, we employ the following MPC strategy:

arg min u(•) Np k=0 ∥x(k)∥ 2 Q + ∥u(k)∥ 2 R + ∥x(N p )∥ 2 P (21a) s.t x(k + 1) = x(k) + τ [f (x(k)) + g(x)u(k)] , u(k) ∈ U, x(k) ∈ X , x(N p ) ∈ X f , (21b)
where Q, R is the user-defined weightings penalizing the tracking and input, respectively. P and X f denote the weighting for the terminal cost and its terminal set, which are associated with a locally optimal control u opt (x). This standard penalization and constraints guarantee the recursive feasibility and stability of the system [START_REF] Mayne | Constrained model predictive control: Stability and optimality[END_REF]. Then, for a sampled state vector, x j , the admissible control u j is taken as the first component of the minimizing sequence of (21b) with x(k = 0) = x j . This formulation hence completes the proposed procedure to collect admissible pairs.

B. Associated nonlinear control law

Prior to this point, the ingredients for the offline generation of the CLF have been presented. Yet, a further question is how one can impose a desired performance over such a CLF during online implementation. Following the continuous stabilizer in [START_REF] Grammatico | Control-sharing and merging control lyapunov functions[END_REF], consider the controller:

u = arg min u ∥u -u d (x)∥ 2 2 (22a) s.t u ∈ U ∇V (x) [f (x) + g(x)u] ≤ -βV (x). (22b)
Namely, the constrained control ensures that the time derivative of V (x) is negative, while the objective function minimizes the deviation of u from the desired one, u d (x). The above optimization is convex due to the convexity of U in (2) and the linear constraint in (22b). From a practical viewpoint, the role of V (x) is to ensure stability, while choosing u d (x) is also important for performance. This control is typically selected as the locally optimal one. Remark 3: Suppose that there exists a DOA D κ = {x : V (x) ≤ κ} deduced from the generated Lyapunov function V (x). It follows directly from the differential inequality (3) that the control [START_REF] Vitale | Lp metrics for compact, convex sets[END_REF] ensures the exponential convergence with the decay constant β inside D κ . Moreover, note that we can always ensure that such a κ exists if we insert among the basis function a locally optimal quadratic function associated with the standard LQ control.

Next, to demonstrate and highlight the applicability of the procedure, let us proceed with a case study of the quadcopter position control problem (a.k.a. outer loop or high level control) in the next section.

IV. APPLICATION ON QUADCOPTER POSITION CONTROL

In this section, via a feedback linearization law, two different representations of the quadcopter's constrained position control problem are presented. Then, the constructions of the CLF for each representation will follow.

A. Problem formulation for quadcopter position control

Let us provide the outer-loop model for a miniature quadcopter as commonly used in the literature [START_REF] Kim | Accurate modeling and robust hovering control for a quad-rotor vtol aircraft[END_REF]:

ξ = T R E -ge 3 , (23)
where ξ = [x, y, z] ⊤ ∈ R 3 collects the drone's position. In the Euler angles representation, ϕ, θ, ψ respectively denote the roll, pitch and yaw angles. R E denotes the rotation matrix from the body to the global frame, which is computed as:

R E =   cos ϕ sin θ cos ψ + sin ϕ sin ψ cos ϕ sin θ sin ψ -sin ϕ cos ψ cos ϕ cos θ   . g is the gravitational
acceleration and e 3 = [0, 0, 1] ⊤ is the standard basis vector pointing along the altitude z axis. T denotes the normalized thrust provided by the four propellers. In a quasi-stationary flight settings, we assume that the aerodynamic forces are always in the direction of the z axis. Centrifugal force and disturbance (e.g., drag) are neglected. As the outer loop (see Fig. 9), the control objective for system [START_REF] Blanchini | Suboptimal receding horizon control for continuous-time systems[END_REF] is to design the position-stabilizing ϕ, θ angles and thrust T . These three signals are regarded as the input for this control layer, which then will be actuated by the low-level controller with a proper feedback scheme and rotor configuration [START_REF] Nguyen | Flat trajectory design and tracking with saturation guarantees: a nano-drone application[END_REF]. In this setup, the input constraint set for [T, ϕ, θ] ⊤ is described as:

C = {0 ≤ T ≤ T max , |ϕ| ≤ ϕ max , |θ| ≤ θ max } (24)
where T max > g, (ϕ max , θ max ) ∈ (0; π/2) 2 are the inputs' constant bounds. It is well known that system ( 23) is feedback linearizable. Indeed, consider the input transformation:

T = u 2 x + u 2 y + (u z + g) 2 , (25a)
ϕ = arcsin ((u x sin ψ -u y cos ψ)/T ), (25b) θ = arctan ((u x cos ψ + u y sin ψ)/(u z + g)),

(25c) then, under the condition u z ≥ -g, the new dynamics can be exactly linearized as a system of three double integrators:

ṗq = A 2D p q + B 2D u q , q ∈ {x, y, z}, (26)
where

p q = [q, q] ⊤ ; A 2D = [ 0 1 0 0 ], B 2D = [0 1] ⊤ .
Besides, the image of the constraint ( 24) via ( 25) can be characterized as a convex set of u pos ≜ [u x , u y , u z ] ⊤ ∈ R 3 as in [START_REF] Do | Indoor experimental validation of mpc-based trajectory tracking for a quadcopter via a flat mapping approach[END_REF]:

u pos ∈ U pos ≜ u pos : ∥u pos + ge 3 ∥ 2 2 ≤ T 2 max ; u 2 x + u 2 y ≤ (u z + g) tan ϵ max . (27)
with ϵ max ≜ min(ϕ max , θ max ). The set U pos is convex since it is the intersection of two convex sets: a ball of radius T max centered at -ge 3 and a convex cone defined with ϵ max and g. A detailed characterization of U pos as in ( 27) can be found in [START_REF] Do | Indoor experimental validation of mpc-based trajectory tracking for a quadcopter via a flat mapping approach[END_REF]. With the dynamics [START_REF] Grammatico | Control-sharing and merging control lyapunov functions[END_REF], one can consider this linearized system as either three independent double integrators in R 2 tangled by the constraint [START_REF] Kim | Accurate modeling and robust hovering control for a quad-rotor vtol aircraft[END_REF] or one system in R 6 and the input constraint u pos ∈ U pos . Next, we investigate the dynamics from both points of view with the proposed procedure. To this aim, we consider two common inner approximations of U pos as in the following. Case 1: A box-type subset: Since U pos includes 0 as an interior point, there always exist ūi > 0 such that:

U b ≜ {u pos ∈ R 3 : |u q | ≤ ūq , q ∈ {x, y, z}} ⊂ U pos . (28)
In this setting, system (26) is decoupled as three linear subsystems associated with three input constraints |u q | ≤ ūq . The maximum volume U b can be found by solving a convex optimization problem, maximizing its volume while constraining all the vertices inside U pos .

Case 2: A polytopic subset: In a less conservative manner, by parameterizing the boundary of U pos , we can obtain a polytopic approximation as:

U p ≜ co [0, 0, -g] ⊤ , [R * cos γ, R * sin γ, u * z ] ⊤ [r cos γ, r sin γ, T 2 max -r 2 -g] ⊤ (29)
with R * = T max sin ϵ max , u * z = T max cos ϵ max -g and for some large integers S 1 , S 2 > 2, γ ∈ N (S 1 , 0, 2π), r ∈ N (S 2 , 0, R * ). For this scenario, we arrange system [START_REF] Grammatico | Control-sharing and merging control lyapunov functions[END_REF] as:

ζ = A 6D ζ + B 6D u pos , (30)
where ζ ≜ p x p y p z ∈ R 6 and u pos denote the state and the input vector as in [START_REF] Kim | Accurate modeling and robust hovering control for a quad-rotor vtol aircraft[END_REF], respectively. u pos ∈ U p as in [START_REF] Do | Indoor experimental validation of mpc-based trajectory tracking for a quadcopter via a flat mapping approach[END_REF],

A 6D = diag(A 2D , A 2D , A 2D ) and B 6D = diag(B 2D , B 2D , B 2D ).
in [START_REF] Kim | Accurate modeling and robust hovering control for a quad-rotor vtol aircraft[END_REF] in ( 28) in ( 29) Fig. 2. The inner approximations of Upos as in [START_REF] Kim | Accurate modeling and robust hovering control for a quad-rotor vtol aircraft[END_REF].

TABLE I NUMERICAL PARAMETERS

Symbols

Values Tmax as in [START_REF] Kim | Accurate modeling and robust hovering control for a quad-rotor vtol aircraft[END_REF] 1.45g ≈ 14.22 m/s 2 θmax = ϕmax = ϵmax as in ( 24), [START_REF] Kim | Accurate modeling and robust hovering control for a quad-rotor vtol aircraft[END_REF] 0.1745 rad ūx, ūy, ūz as in [START_REF] Nguyen | Flat trajectory design and tracking with saturation guarantees: a nano-drone application[END_REF] 0.8154, 0.8154, 3.27 S 1 , S 2 as in [START_REF] Do | Indoor experimental validation of mpc-based trajectory tracking for a quadcopter via a flat mapping approach[END_REF] 7, 3

An illustration and numerical parameters of the abovementioned sets are depicted in Fig. 2 and TABLE I, respectively.

Remark 4: It is undoubted that, Procedure 1 can be immediately applied for the system (30) in Case 2. However, we believe that exploiting the decoupled representation in Case 1, i.e. examining each double integrator at a time, provides not only illustrative understanding but also fundamental insights for the choice of the basis functions V i (•) as in [START_REF] Freeman | Inverse optimality in robust stabilization[END_REF]. Thus, in the sequel, we proceed by investigating first the decoupled representation in Case 1 (or ( 26)), then the procedure will be extended to Case 2 with the concatenated dynamics [START_REF] Blanchini | A new class of universal Lyapunov functions for the control of uncertain linear systems[END_REF]. □

B. CLF construction

In this subsection, the formulation of the CLFs for ( 26) and [START_REF] Blanchini | A new class of universal Lyapunov functions for the control of uncertain linear systems[END_REF] subject to [START_REF] Nguyen | Flat trajectory design and tracking with saturation guarantees: a nano-drone application[END_REF] and [START_REF] Do | Indoor experimental validation of mpc-based trajectory tracking for a quadcopter via a flat mapping approach[END_REF], respectively, will be discussed. With these functions come two ways of controlling the drone, as the following. In the former, three CLFs will be generated for each axis q ∈ {x, y, z}. These functions then will be separately combined with the setup [START_REF] Vitale | Lp metrics for compact, convex sets[END_REF] to compute the controller u q . In the latter, one CLF will be generated for the concatenated dynamics [START_REF] Blanchini | A new class of universal Lyapunov functions for the control of uncertain linear systems[END_REF] with the polytopic constraint [START_REF] Do | Indoor experimental validation of mpc-based trajectory tracking for a quadcopter via a flat mapping approach[END_REF]. Then, the controller u pos = [u x , u y , u z ] ⊤ will be computed similarly with [START_REF] Vitale | Lp metrics for compact, convex sets[END_REF] and the generated function. In both cases, the controller [u x , u y , u z ] ⊤ will be transformed back to the inputs T, ϕ, θ by means of [START_REF] Mayne | Constrained model predictive control: Stability and optimality[END_REF], then applied to the vehicle.

For Case 1: The idea is to apply Procedure 1 to three subsystems in [START_REF] Grammatico | Control-sharing and merging control lyapunov functions[END_REF] with |u q | ≤ ūq , q ∈ {x, y, z} and to create one CLF for each double integrator. For simplicity, we use V q (•) to denote the generated CLF for the q axis (e.g., V x (p), p ∈ R 2 , is the CLF for the first double integrator in [START_REF] Grammatico | Control-sharing and merging control lyapunov functions[END_REF], or q = x). Then, the same procedure is carried out for all three subsystems as follows.

First, the general admissible pair generator ( 21) will be adapted with the weighting as in The domain of interest for each double integrator is chosen as X q = {|p q | ≤ 1.5}, q ∈ {x, y, z}.

The next step resides on how to choose a collection basis functions V i (•) as in [START_REF] Leyva | Stabilization of affine systems with polytopic control value sets[END_REF]. On one hand, encoding the optimality, at least locally near the origin, is important. For this reason, we select the first basis as the quadratic function:

V 1 (p) = ∥p∥ 2 P lqr 2D , p ∈ R 2 . ( 31
)
where P lqr 2D denotes the result of the continuous Riccati equation, with the state and input penalty abusively taken from Q 2D and R . On the other hand, it is necessary that the basis functions are sufficiently rich so that ( 6) is feasible. At the same time, it was proposed that a suitable Lyapunov function for constrained stabilization is the 2p-norm ∥Gp∥ 2p 2p where for p ∈ R n , G ∈ R n×r is a full row rank matrix [START_REF] Blanchini | A new class of universal Lyapunov functions for the control of uncertain linear systems[END_REF]. Thus, the remaining members of the collection of V i (p) are particularly populated with the following polynomials:

V 2≤i≤100 (p) ∈ ∥ηF p∥ 4 4 , η ∈ ls(33, 0.4, 2.0)

F ∈ 1 0 , 0 1 , 1 1 . ( 32
)
Remark 5: In (32), η plays a role of increasingly enriching the basis with different value scale, while F is responsible for the "orientation" of the surfaces V i (•). The 2D case of the first double integrator is depicted in Fig. 3. The feasibility of ( 6) can be managed by meticulously collecting more value of η and vector F . □ Then, by partitioning the domain X q by a 25 × 25 grid and collecting the control u q given by ( 21) at those points (See Fig. 6), we achieved the three following CLFs from [START_REF] Freeman | Inverse optimality in robust stabilization[END_REF] with the basis ( 31) and (32):

• For q ∈ {x, y} in (26):

V q (p q ) = 0.016∥p q ∥ 2

P lqr 2D + 0.2293(2[1 0]p q ) 4 +0.5406(2[0 1]p q ) 4 + 0.2141(2[1 1]p q ) 4 . (33)
• For the z axis in ( 26):

V z (p z ) = 0.0119∥p z ∥ 2 P lqr 2D + 0.3862(2[1 0]p z ) 4 +0.6019(2[1 1]p z ) 4 (34)
V x (p x )
as in [START_REF] Förster | System identification of the crazyflie 2.0 nano quadrocopter[END_REF] Fig. 4. Generated Lyapunov functions and its DOA for the first subsystem of (26) (q = x).

-1.5 -1 -0.5 0 0.5 1 1.5

x (m) time (s)

-1 0 1 input ux V x (p x ) + βV x (p x ) V x (p x )
Fig. 5. Set-point tracking simulation with the first subsystem of (26) (q = x). The green-dashed line represents the DOA provided by the MPC point generator [START_REF] Hosseinzadeh | Robust to early termination model predictive control[END_REF]. The red solid curve depicts the sublevel set DOA given by the generated Lyapunov function.

The computation of α i in ( 6) is carried out with MATLAB 2021b and its linprog solver, while the pair generator [START_REF] Hosseinzadeh | Robust to early termination model predictive control[END_REF] and the online controller ( 22) is computed with quadprog solver for the simulation. The illustration for ( 33) is in Fig. 4. The DOA can be estimated by verifying the state-space with the Fig. 6. Comparison between the control sampled from (21) (blue points) and the control ( 22) (pink surface) for Case 1, q = x. condition [START_REF] Prajna | Nonlinear control synthesis by sum of squares optimization: A lyapunov-based approach[END_REF], which, in this case, is given as the sublevel set V x (p x ) ≤ 23.04. The positive invariance of the DOA can also be confirmed via the simulated trajectories and vector field depicted in Fig. 5. Therein, it can be seen again that there always exists |u x | ≤ ūx rendering V x (p x ) decreasing and V x (p x ) + βV x (p x ) ≤ 0. A comparison between the control policy of the MPC ( 21) and ( 22) is given in Fig. 6 .

For Case 2: Similarly, we use the setting [START_REF] Hosseinzadeh | Robust to early termination model predictive control[END_REF] as the point generator with the concatenated state ζ as in [START_REF] Blanchini | A new class of universal Lyapunov functions for the control of uncertain linear systems[END_REF] and the parameters computed as

Q 6D = diag(Q 2D , Q 2D , Q 2D ), R 6D = diag(R 2D , R 2D , R 2D
), P 6D = diag(P 2D , P 2D , P 2D ) and consequently, K 6D = diag(K 2D , K 2D , K 2D ). Parameters' details are given in TABLE II.

Different from the previous case, for the admissible point collection, random initial conditions in {ζ ∈ R 6 : |ζ| ≤ 1.5}. The toolbox given in [START_REF] Benham | Uniform distribution over a convex polytope[END_REF] is employed to provide a randomly distributed initial conditions. Then, for each feasible initial condition, we sample along its trajectory to collect statecontrol pairs since the recursive feasibility of (21) ensures the existence of admissible pairs propagated from that initial state. The projections of those states are provided in Fig. 7. Then the new set of basis functions for this setup is reselected as follows: The first basis function remains as the LQR-based quadratic function in the same choice from (31):

V 1 (ζ) = ∥ζ∥ 2 P lqr , ζ ∈ R 6 . ( 35
)
Moreover, taking advantage of the successfully generated CLFs for each subspace in ( 33) and ( 34), we adopt their basis by augmenting the vector F in (32), for q ∈ {x, y, z}, with the extended value assigned as 0. For example, for q = x, [1 1] becomes [1 1 0 0 0 0]. For the tuning of ℓ i , as suggested before, it is essential to include the LQR-based quadratic Lyapunov function in the parameterized CLF in order to attain local optimality. Therefore, in both cases, the weighting ℓ 1 associated with the function is progressively decreased so that its corresponding coefficient computed from ( 6) is non-zero.

With the pairs sampled and the choice of basis functions given, the LP (6) provides the CLF as: (36) with V (ζ) denoting the generated Lyapunov function for the six-dimensional system [START_REF] Blanchini | A new class of universal Lyapunov functions for the control of uncertain linear systems[END_REF]. The computational details of the procedure are provided in TABLE III. Therein, as expected, it is noticeable that the computation time for the LP ( 6) is relatively small in both cases, confirming the low complexity of the technique.

V (ζ) = 0.
To this end, the CLFs were generated for the control problem and examined via simulation tests. In the next section, we verify the constructed CLFs via experimental validation, showing their applicability in practice.

V. EXPERIMENTAL VALIDATION

Herein, the results from the previous section are validated via experimental tests over nano-drone platforms.

Motion capture cameras Active markers Fig. 8. The Crazyflie's motion restored with a system of cameras capturing the infrared lights from the active markers.

A. Experiment setup

In this section, the open-source quadcopter, Crazyflie 2.1, is used, since it allows customizing the position control within its hierarchical structure (See Fig. 9). We further assume that the integrated inner loop (the green block) is capable of sustaining the rotational dynamics in a properly high bandwidth compared to that of the outer loop, thereby ensuring that the roll and pitch angle commands from the outer loop are tracked with sufficient speed. A detailed description of this layer can be found in [START_REF] Huang | Closed-loop identification and realtime control of a micro quadcopter[END_REF]. For the outer loop, we substitute our position control, which is designed for system (26) (or equivalently [START_REF] Blanchini | A new class of universal Lyapunov functions for the control of uncertain linear systems[END_REF]) into the blue block to track the reference position ξ ref .

The control action u pos then transformed back to the inputs T, ϕ and θ by means of [START_REF] Mayne | Constrained model predictive control: Stability and optimality[END_REF]. Note that the normalized thrust T needs to be converted into the corresponding value in the digital scale of the nanodrone, ranging from 0 to 65535. The conversion is done via an interpolated curve identified experimentally (see Chapter 3,[START_REF] Förster | System identification of the crazyflie 2.0 nano quadrocopter[END_REF]). These signals are then packed with the value of the desired yaw rate, ψref ≜ 0. These calculations are carried out in a station computer and and sent to the onboard inner loop via a USB radio dongle. The sampling time is fixed at t s = 100 ms. For all the optimization problems employed in the experiments, the IPOPT solver is used within the Casadi optimization framework in Python 3.8.8, Intel Core, i5-10300H CPU @ 2.50GHz and 16GB RAM. The description of the test cases is given as follows.

On one hand, we use the controller [START_REF] Vitale | Lp metrics for compact, convex sets[END_REF] with the two setups Case 1 and Case 2 in the previous section with the same tuning given in TABLE II

N * p k=0 ∥ζ(k) -ζ ref (k)∥ 2 Q 6D + ∥u pos (k)∥ 2 R 6D + ∥ζ(N * p ) -ζ ref (N * p )∥ 2 P 6D (37a)      ζ(k + 1) = (A 6D + t s I)ζ(k) + t s B 6D u pos (k), u pos (k) ∈ U pos as in (29), ζ(k) -ζ ref (k) ∈ Z, ζ(N * p ) -ζ ref (N * p ) ∈ Z f , (37b)
where

A 6D , B 6D , P 6D , Q 6D , R 6D are given in Section IV-A, ζ ref (k) denotes the reference for ζ(k) to follow. Z = {ζ : ∥ζ∥ ∞ ≤ 1.
5} is the state constraint corresponding to that of the admissible pair generator in Case 1 and 2. Z f is the MPI set (as the terminal constraint) associated with the full-state gain K 6D given therein. For our setup, the largest prediction horizon N * p for (37) is 20, so that the online computation time does not exceed the sampling time t s .

On the other hand, for the position reference, the examined controllers will first track a set-point ξ ref = [0.6, 0.6, 0.8] ⊤ (m), noted as Scenario 1 (Sce. 1). Then, to further emphasize the computational advantage, we apply the same control laws for three Crazyflies in a centralized manner. The control then will be sent to each drone via the Swarm Python interface [START_REF]Step-by-step: Swarm interface[END_REF]. With this Scenario 2 (Sce. 2), a position swapping reference was given to the three drones, which is depicted Fig. 12.

B. Experimental results and discussions

The result of the first scenario (Sce. 1) is shown in Fig. 10. It can be seen that, the proposed schemes' tracking performances (Case 1 and 2) are commensurate with the MPC while respecting all the input constraints in [START_REF] Rockafellar | Convex analysis[END_REF]. However, it is noted from TABLE IV and Fig. 14 that the computation time (CT) of the proposed laws is 5 to 7 times lower than that of the MPC (37). This advantage directly comes from the simple QP formulation of (22) compared to (37) with regard to the complexity of both the cost and constraints. In Sce. 2, the proposed schemes are also successfully validated. It is noticeable that in this scenario, a centralized MPC for three quadcopters could not be performed since the required calculation time exceeds the sampling time t s . Meanwhile, even with three quadcopters, the tracking was possible for both Case 1 and 2 with the Avg. CT around 20ms. Fig. 11 depicts the value of the CLF in both cases and their monotonicity. This shows that the discussed technique generates effective CLF and is applicable in real-life implementation.

VI. CONCLUSION

In this work, an offline synthesis procedure for CLF (Control Lyapunov Function) was proposed via linear programming. The idea relied on the stabilizing MPC (Model Predictive Control)-based software to generate admissible control in the state-space and parameterizing a collection of candidate CLFs. The generated CLF's advantage was then validated via quadcopter stabilization experiments. The result showed that the offline calculation of the CLF can be implemented via a

  terminal cost adopts the weighting P 2D from the discrete Riccati equation for the EAS of (26) (τ = 0.1s) with the penalty weighting Q 2D , R 2D , while the local controller uses the standard infinite-horizon unconstrained control. The terminal region X f is the maximal positive invariant (MPI) set computed with such a local control and input constraint |u q | ≤ ūq for each subsystem. Numerical details are given as Q 2D = [ 50 0 0 5 ], R 2D = 5, P 2D = [ 479.6118 181.0469 181.0469 155.5598 ], K 2D = -[2.7617 2.6491].

Fig. 3 .

 3 Fig. 3. The surface of V i (•) with different η and F (Case 1).

Fig. 7 .

 7 Fig. 7. Admissible pairs (red) collected for Case 2 (projected to each positionvelocity axis).

0244∥ζ∥ 2 P

 2 lqr 6D + 0.4508([0 1 0 0 0 0]ζ) 4 + 0.0674([1 1 0 0 0 0]ζ) 4 + 0.1631([0 0 1 0 0 0]ζ) 4 + 0.1841([0 0 1 1 0 0]ζ) 4 + 0.1103([0 0 0 0 0 1]ζ) 4 ,

  . The desired control u d (•) is chosen as u d (p q ) = -[k p , k d ]p q and u d (ζ) = -diag([k p , k d ], [k p , k d ], [k p , k d ])ζ for Case 1 and 2, respectively, and k p = 1, k d = 1.45. Then, for comparison, we stabilize the tracking error dynamics by employing the following optimization problem in a receding horizon fashion.

  min

Fig. 9 .

 9 Fig. 9. Hierarchical control architecture of the Crazyflie quadcopter.

Fig. 10 .

 10 Fig. 10. Left: Set-point reference tracking experimented with the three scenarios. Right: The corresponding real input signals and their amplitude bound (dashed-dotted line).

Fig. 11 .Fig. 12 .

 1112 Fig. 11. Lyapunov function for Case 1 and 2 controllers in Scenario 1.

Fig. 13 .

 13 Fig. 13. The outer loop's input signals of the three quadcopters in Sce. 2 with the Case 2 controller.

Fig. 14 .

 14 Fig. 14. Computation time in Scenario 1.

TABLE II PARAMETERS

 II FOR POINT GENERATOR[START_REF] Hosseinzadeh | Robust to early termination model predictive control[END_REF] WITH τ = 0.1s, AND THE LP (6), q ∈ {x, y, z}.

	Representation	Variables x, u	Q	R	P	Constraints X , U for (21)	Np	β	ℓ i
	Case 1	p x , ux p y , uy p z , uz	Q 2D	R 2D	P 2D	|p q | ≤ 1.5, |uq| ≤ ūq	50	0.1275 0.1275 2.65	ℓ 1 = 0.05 ℓ i̸ =1 = 0.1
	Case 2	ζ, upos	Q 6D	R 6D	P 6D	|ζ| ≤ 1.5, upos ∈ Upos	50	0.15	ℓ 1 = 0.2 ℓ i̸ =1 = 1

TABLE III SPECIFICATIONS

 III OF THE LP (6) FOR THE TWO CASES

			Case 1		Case 2
		V x (p x )	V y (p y )	V z (p z )	V (ζ)
	No. basis functions N	100	100	100	10
	No. admissible pairs M	625	625	625	6330
	Computation time of (6)	0.03 (s)	0.03 (s)	0.05 (s)	0.11 (s)

  linear program, thereby allowing CLF construction for systems with nontrivial dimensions. With the experimental results, the generated non-quadratic CLF, can be employed to reduce the online computational demand while maintaining stability. The future work will investigate how to efficiently choose the sampled data from MPC, hence, providing theoretical guarantees over a finite set of constraints, and allowing scalability.

			TABLE IV
			EXPERIMENT RESULTS
		Controller	RMS tracking error (cm)	Avg. CT
		Case 1		14.11		7.02 ms
	Sce. 1	Case 2		11.39		10.74 ms
		MPC		11.48		54.70 ms
			Drone 1	Drone 2	Drone 3
	Sce. 2	Case 1 Case 2	24.09 19.14	22.20 18.57	21.33 18.71	24.09 ms 20.51 ms
		MPC	×	×	×	×