Ségolène Martin

Malik Boudiaf

Éts Montreal

Emilie Chouzenoux

Jean-Christophe Pesquet

Ismail Ben Ayed

Towards Practical Few-Shot Query Sets: Transductive Minimum Description Length Inference

Standard few-shot benchmarks are often built upon simplifying assumptions on the query sets, which may not always hold in practice. In particular, for each task at testing time, the classes effectively present in the unlabeled query set are known a priori, and correspond exactly to the set of classes represented in the labeled support set. We relax these assumptions and extend current benchmarks, so that the query-set classes of a given task are unknown, but just belong to a much larger set of possible classes. Our setting could be viewed as an instance of the challenging yet practical problem of extremely imbalanced K-way classification, K being much larger than the values typically used in standard benchmarks, and with potentially irrelevant supervision from the support set. Expectedly, our setting incurs drops in the performances of state-of-the-art methods. Motivated by these observations, we introduce a PrimAl Dual Minimum Description LEngth (PADDLE) formulation, which balances data-fitting accuracy and model complexity for a given few-shot task, under supervision constraints from the support set. Our constrained MDL-like objective promotes competition among a large set of possible classes, preserving only effective classes that befit better the data of a few-shot task. It is hyperparameter free, and could be applied on top of any base-class training. Furthermore, we derive a fast block coordinate descent algorithm for optimizing our objective, with convergence guarantee, and a linear computational complexity at each iteration. Comprehensive experiments over the standard few-shot datasets and the more realistic and challenging i-Nat dataset show highly competitive performances of our method, more so when the numbers of possible classes in the tasks increase.

Introduction

The performance of deep learning models is often seriously affected when tackling new tasks with limited supervision, i.e., classes that were unobserved during training and for which we have only a handful of labeled examples. Few-shot learning [START_REF] Snell | Prototypical networks for few-shot learning[END_REF][START_REF] Finn | Model-agnostic meta-learning for fast adaptation of deep networks[END_REF][START_REF] Vinyals | Matching networks for one shot learning[END_REF] focuses on this generalization challenge, which occurs in a breadth of applications. In standard few-shot settings, a deep network is first trained on a large-scale dataset including labeled instances sampled from an initial set of classes, commonly referred to as the base classes. Subsequently, for novel classes, unobserved during the base training, supervision is restricted to a limited number of labeled instances per class. During the test phase, few-shot methods are evaluated over individual tasks, each including a small batch of unlabeled test samples (the query set) and a few labeled instances per novel class (the support set).

The transduction approach has becomes increasingly popular in few-shot learning, and a large body of recent methods focused on this setting, including, for instance, those based on graph regularization [START_REF] Ziko | Laplacian regularized few-shot learning[END_REF][START_REF] Yanbin | Learning to propagate labels: Transductive propagation network for few-shot learning[END_REF], optimal transport [START_REF] Lazarou | Iterative label cleaning for transductive and semi-supervised few-shot learning[END_REF][START_REF] Hu | Leveraging the feature distribution in transfer-based few-shot learning[END_REF], feature transformations [START_REF] Liu | Prototype rectification for few-shot learning[END_REF][START_REF] Cui | Parameterless transductive feature re-representation for few-shot learning[END_REF], information maximization [START_REF] Veilleux | Realistic evaluation of transductive few-shot learning[END_REF][START_REF] Boudiaf | Transductive information maximization for few-shot learning[END_REF][START_REF] Boudiaf | Few-shot segmentation without meta-learning: A good transductive inference is all you need?[END_REF] and transductive batch normalization [START_REF] Bronskill | Tasknorm: Rethinking batch normalization for meta-learning[END_REF][START_REF] Finn | Model-agnostic meta-learning for fast adaptation of deep networks[END_REF], among other works [START_REF] Dhillon | A baseline for few-shot image classification[END_REF][START_REF] Wang | Instance credibility inference for few-shot learning[END_REF][START_REF] Cui | Parameterless transductive feature re-representation for few-shot learning[END_REF][START_REF] Qi | Transductive few-shot classification on the oblique manifold[END_REF][START_REF] Shen | Re-ranking for image retrieval and transductive few-shot classification[END_REF]. Transductive few-shot classifiers make joint predictions for the batch of query samples of each few-shot task, independently of the other tasks. Unlike inductive inference, in which prediction is made for one testing sample at a time, transduction accounts for the statistics of the query set of a task 4 , typically yielding substantial improvements in performance. On standard benchmarks, the gap in classification accuracy between transductive and inductive few-shot methods may reach 10%; see [START_REF] Boudiaf | Transductive information maximization for few-shot learning[END_REF], for example. This connects with a well-known fact in classical literature on transductive inference [START_REF] Vapnik | An overview of statistical learning theory[END_REF][START_REF] Joachims | Transductive inference for text classification using support vector machines[END_REF][START_REF] Dengyong | Learning with local and global consistency[END_REF], which prescribes transduction as an effective way to mitigate the difficulty inherent to learning from limited labels. It is worth mentioning that transductive methods inherently depend on the statistical properties of the query sets. For instance, the recent studies in [START_REF] Lichtenstein | TAFSSL: Task-adaptive feature sub-space learning for few-shot classification[END_REF][START_REF] Veilleux | Realistic evaluation of transductive few-shot learning[END_REF] showed that variations in the class balance within the query set may affect the performances of transductive methods.

Transductive few-shot classification occurs in a variety of practical scenarios, in which we naturally have access to a batch of unlabeled samples at test time. In commonly used few-shot benchmarks, the query set of each task is small (less than 100 examples), and is sampled from a limited number of classes (typically 5). Those choices are relevant in practice: During test time, one typically has access to small unlabeled batches of potentially correlated (non-i.i.d.) samples, e.g., smart-device photos taken at a given time, video-stream sequences, or in pixel prediction tasks such as semantic segmentation [START_REF] Boudiaf | Few-shot segmentation without meta-learning: A good transductive inference is all you need?[END_REF], where only a handful of classes appear in the testing batch. However, the standard few-shot benchmarks are built upon further assumptions that may not always hold in practice: the few classes effectively present in the unlabeled query set are assumed both to be known beforehand and to match exactly the set of classes of the labeled support set. We relax these assumptions, and extend current benchmarks so that the query-set classes are unknown and do not match exactly the support-set classes, but just belong to a much larger set of possible classes. Specifically, we allow the total number of classes represented in the support set to be higher than typical values, while keeping the number of classes that are effectively present in the query set to be much smaller.

Our challenging yet practical setting raises several difficulties for state-of-the-art transductive fewshot classifiers: (i) it is an instance of the problem of highly imbalanced classification; (ii) the labeled support set includes "distraction" classes that may not actually be present in the test query samples; (iii) we consider K-way classification tasks with K much larger than the typical values used in the current benchmarks. We evaluated 7 of the best-performing state-of-the-art transductive few-shot methods and, expectedly, observed drops in their performance in this challenging setting.

Motivated by these experimental observations, we introduce a Minimum Description Length (MDL) inference, which balances data-fitting accuracy and model complexity for a given few-shot task, subject to supervision constraints from the support set. The model-complexity term can be viewed as a continuous relaxation of a discrete label cost, which penalizes the number of non-empty clusters in the solution, fitting the data of a given task with as few unique labels as necessary. It encourages competition among a large set of possible classes, preserving only those that fit better the task. Our formulation is hyper-parameter free, and could be applied on top of any base-class training. Furthermore, we derive a fast primal-dual block coordinate descent algorithm for optimizing our objective, with convergence guarantee, and a linear computational complexity at each iteration thanks to closed-form updates of the variables. We report comprehensive experiments and ablation studies on mini-Imagenet, tiered-Imagenet, and the more realistic and challenging iNat dataset [START_REF] Wertheimer | Few-shot learning with localization in realistic settings[END_REF] for fine-grained classification, with 908 classes and 227 ways at test-time. Our method yields competitive performances in comparison to the state-of-the-art, with gaps increasing substantially with the numbers of possible classes in the tasks.

Few-shot setting and task generation

Base training Let D base = {x n , y n } |Dbase| n=1 denotes the base dataset, where each x n ∈ X base is a sample from some input space X , y n ∈ {0, 1} |Ybase| the associated ground-truth label, and Y base the set of base classes. Base training learns a feature extractor f φ : X → Z, with parameters φ and Z a lower-dimensional space. For this stage, an abundant few-shot literature adopts episodic training, which views D base as a series of tasks (or episodes) so as to simulate testing time. Then, a meta-learner is devised to produce the predictions. However, it has been widely established over recent years that a basic training, followed by transfer-learning strategies, outperforms most meta-learning methods [START_REF] Chen | A closer look at few-shot classification[END_REF][START_REF] Wang | Simpleshot: Revisiting nearest-neighbor classification for few-shot learning[END_REF][START_REF] Tian | Rethinking few-shot image classification: a good embedding is all you need?[END_REF][START_REF] Ziko | Laplacian regularized few-shot learning[END_REF][START_REF] Boudiaf | Transductive information maximization for few-shot learning[END_REF]. Hence, we adopt a standard cross-entropy training in this work. Task generation For a given task, let K denote the total number of possible classes that are represented in the labeled support set S, and K eff the number of classes that appear effectively in unlabeled query set Q (i.e. classes represented by at least one sample). The standard task generation protocol assumes that the set of effective classes in Q matches exactly the set of classes in S, i.e. K = K eff |Y test |, which amounts to knowing exactly the few classes that appear in the test samples. Then, a fixed number of instances per class are sampled for each query set, forcing class balance. We relax these assumptions, and propose a setting where the set of effective classes in Q is not known exactly and belongs to a much larger set of possible classes. First, we allow the total number of possible classes K to be higher than typical values, with K = |Y test |, while keeping the number of effective classes in the query set to be typically much smaller: K eff |Y test | = K. For instance, in our experiments, K can go up to 227. This, as expected, results in K-way problems that are more challenging than the standard 5-way tasks used in the few-shot literature. Secondly, once the effective classes are fixed, and given a budget of images, we sample the query set under the data joint distribution (i.e. uniform draws among all available examples), so that the statistics of the sampled query sets more faithfully reflect the natural distribution of classes. Figure 1 illustrates our framework (base training, inference and task generation).

Proposed few-shot inference formulation

For a given few-shot task, let

Q = {x n } n∈I Q and S = {x n , y n } n∈I S be two subsets of D test such that Q ∩ S = ∅. Let N = |Q| + |S|.
Up to a reordering, we can suppose that I Q = {1, . . . , |Q|} and I S = {|Q| + 1, . . . , N }. We denote z n = f φ (x n ) the feature vector corresponding to the data sample x n , and u n = (u n,k ) 1≤k≤K ∈ {0, 1} K the variable assigning the data to one of the possible classes in {1, . . . , K}, i.e. u n,k = 1 if x n belongs to class k and 0 otherwise. We define the variables giving the class proportions û

= (û k ) 1≤k≤K ∈ ∆ K as ûk = 1 |Q| |Q| n=1 u n,k ∀k ∈ {1, . . . , K}, (1)
where ∆ K is the unit simplex of R K .

We propose to cast transductive few-shot inference as the minimization of an objective balancing data-fitting accuracy and partition complexity, subject to supervision constraints (known labels) from the support set S. Our approach estimates jointly U = (u n ) 1≤n≤N , which defines a partition of The support set includes "distraction" classes that may not actually be present in the query set, e.g. classes "Alligator", "Peacock" and "Butterfly". All possible classes (K classes) are represented in the support set, but only an unknown subset among these K possible classes (K eff effective classes) appear in the query set, with K eff K.

S ∪ Q, and the class prototypes

W = (w k ) 1≤k≤K ∈ (R d ) K through the following problem: minimize U ,W 1 2 K k=1 N n=1 u n,k w k -z n 2 data-fitting accuracy -λ K k=1 ûk ln(û k ) partition complexity , (2)
s.t u n ∈ ∆ K ∀n ∈ {1, . . . , |Q|}, u n = y n ∀n ∈ {|Q| + 1, . . . , N }. (C)
Note that, in the second line of (2), we have relaxed the integer constraints on query assignments to facilitate optimization.

Effect of each term

The purpose of objective ( 2) is to classify the data of a few-shot task with as few unique labels as necessary. The data-fitting term has the same form as the standard K-means objective for clustering, but is constrained with supervision from the support-set labels. This term evaluates, within each class k, the deviation of features from class prototype w k , thereby encouraging consistency of samples of the same class. The partition-complexity term implicitly penalizes the number of effective (non-empty) classes that appear in the solution (K eff ≤ K ), encouraging low cardinality partitions of query set Q. This term is the Shannon entropy of class proportions within Q: it reaches its minimum when all the samples of Q belong to a single class, i.e., ∃j ∈ {1, . . . , K} such that ûj = 1 and all other proportions vanish. It achieves its maximum for perfectly balanced partitions of Q, satisfying ûk = 1/K for all k. In practice, this terms promotes solutions that contain only a handful of effective classes among a larger set of K possible classes.

Connection to Minimum Description Length (MDL)

The objective in (2) could be viewed as a partially-supervised instantiation of the general MDL principle. Originated in information theory, MDL is widely used in statistical model selection [START_REF] Mackay | Information Theory, Inference, and Learning Algorithms[END_REF]. Assume that we want to describe some input data Z = (z n ) 1≤n≤|Q| with a statistical model M , among a family of possible models. MDL prescribes that the best model corresponds to the shortest description of the data, according to some coding scheme underlying the model. It balances data-fitting accuracy and model complexity by minimizing L(Z|M ) + λL(M ) w.r.t M . L(Z|M ) measures the code length of the prediction of Z made by M , while L(M ) encourages simpler models (Occam's razor [START_REF] Mackay | Information Theory, Inference, and Learning Algorithms[END_REF]), e.g., models with less parameters. For instance, in unsupervised clustering, it is common to minimize a discrete label cost as a measure of model complexity [START_REF] Delong | Fast approximate energy minimization with label costs[END_REF][START_REF] Yuan | Tv-based multi-label image segmentation with label cost prior[END_REF][START_REF] Zhu | Region competition: unifying snakes, region growing, and bayes/mdl for multiband image segmentation[END_REF]. Such a label cost is the number of effective (non-empty) clusters in the solution, and is used in conjunction with a log-likelihood term for data-fitting:

L(M ) = K k=1 1 u k =0 (3)
L(Z|M ) = - K k=1 |Q| n=1 u n,k ln Pr(z n |k; w k ) (4)
In this expression, M = {(w k ) 1≤k≤K , U } with w k being the parameters of some probability distribution Pr(z n |k; w k ) describing the samples in cluster k, while U contains assignments as denoted above.

The discrete model complexity in ( 3) is commonly used in MDL-based clustering [START_REF] Delong | Fast approximate energy minimization with label costs[END_REF][START_REF] Yuan | Tv-based multi-label image segmentation with label cost prior[END_REF][START_REF] Zhu | Region competition: unifying snakes, region growing, and bayes/mdl for multiband image segmentation[END_REF], despite the ensuing optimization difficulty. As this measure is a discrete count of non-empty clusters, it does not accommodate fast optimization techniques. Typically, it is either handled via cluster merging heuristics (starting from an initial set of clusters) [START_REF] Zhu | Region competition: unifying snakes, region growing, and bayes/mdl for multiband image segmentation[END_REF] or via combinatorial move-making algorithms [START_REF] Delong | Fast approximate energy minimization with label costs[END_REF], both of which are computationally intensive, more so when dealing with large sets of classes.

Our partition-complexity term in (2) can be viewed as a continuous relaxation of the discrete label count in (3) (in the supplemental material, we provide a graphical illustration in the case when K = 2).

While penalizing similarly the number of effective classes, it has several advantages over the label count in (3), despite the surprising fact that it is not common in the MDL-based clustering literature, to our best knowledge. First, it is a continuous function of assignment variables U . This enables us to derive a fast primal-dual block coordinate descent algorithm (Section 4), with convergence guarantee, and a linear computational complexity at each iteration thanks to closed-form updates of the variables. Second, it has a clear MDL interpretation: it measures the number of bits required to encode the set of classes, given class probabilities Pr(k) = ûk , ∀k ∈ {1, . . . , K}.

Following the Kraft-McMillan theorem [START_REF] Mackay | Information Theory, Inference, and Learning Algorithms[END_REF], any probability distribution Pr(z n |k; w k ) corresponds to some coding scheme for storing the features of cluster k, andln Pr(z n |k; w k ) is the number of bits required to represent any feature using coding scheme Pr(z n |k; w k ). Clearly, our data-fitting term in (2) also fits into this MDL interpretation by assuming that each probability Pr(z n |k; w k ) is a Gaussian distribution with mean w k and covariance fixed to the identity matrix:

Pr(z n |k; w k ) ∝ exp - 1 2 w k -z n 2 ∀k ∈ {1, . . . , K}. (5)
Differences with the existing objectives for transductive few-shot inference Unlike most of the existing transductive few-shot methods, our MDL formulation does not encode strong assumptions on the label statistics of the query set (e.g. class balance or/and a perfect match between the query and support classes). Information maximization methods, such as [START_REF] Veilleux | Realistic evaluation of transductive few-shot learning[END_REF][START_REF] Boudiaf | Transductive information maximization for few-shot learning[END_REF], maximise the confidence of the sample-wise predictions while encouraging class balance. Optimal transport methods [START_REF] Lazarou | Iterative label cleaning for transductive and semi-supervised few-shot learning[END_REF][START_REF] Hu | Leveraging the feature distribution in transfer-based few-shot learning[END_REF] estimate an optimal mapping matrix, which could be viewed as a joint probability distribution over the features and the labels, while imposing a hard class-balance constraint via the Sinkhorn-Knopp algorithm. Both information-maximization and optimal-transport methods have inherent class-balance bias and, as shown in the experiments, undergo a drastic drop in accuracy when the number of possible classes increases (as this corresponds to highly imbalanced problems). Prototype rectification [START_REF] Liu | Prototype rectification for few-shot learning[END_REF] transforms the query features so as to minimize the difference between the overall statistics of the query and support sets. This relies on the assumption that the support and query classes match perfectly. Therefore, it is not adapted to our setting. Inspired from graphical models, the Laplacian regularization in [START_REF] Ziko | Laplacian regularized few-shot learning[END_REF] is a pairwise label correction (rather than a learning) method, which encourages assigning the same label to query samples that are close in the input space; it does not learn a class representation from the query set (the prototypes being fixed as those of the support set).

Identifying λ through an unbiased probabilistic K-means interpretation of (2) Another interesting view of our few-shot inference objective in (2) can be drawn from the probabilistic K-means objective 5 , well-known in the clustering literature [START_REF] Kearns | An information-theoretic analysis of hard and soft assignment methods for clustering[END_REF][START_REF] Tang | Kernel cuts: Kernel and spectral clustering meet regularization[END_REF][START_REF] Boykov | Volumetric bias in segmentation and reconstruction: Secrets and solutions[END_REF]. In fact, probabilistic K-means corresponds to minimizing L(Z|M ) in (4) as a generalization of K-means, which corresponds to the particular Gaussian choice in [START_REF] Yanbin | Learning to propagate labels: Transductive propagation network for few-shot learning[END_REF]. It is a well-known that probabilistic K-means has a strong bias towards balanced partitions [START_REF] Kearns | An information-theoretic analysis of hard and soft assignment methods for clustering[END_REF][START_REF] Boykov | Volumetric bias in segmentation and reconstruction: Secrets and solutions[END_REF]. The objective can be decomposed to reveal a hidden term promoting class balance. Using Bayes rule Pr(k|z n ; w k ) ∝ Pr(z n |k; w k )Pr(k), empirical estimates of marginal class probabilities Pr(k) = ûk , and the choice of Pr(z n |k; w k ) in ( 5), we have:

1 2 K k=1 |Q| n=1 u n,k w k -z n 2 = - K k=1 |Q| n=1 u n,k ln Pr(z n |k; w k ), c = - K k=1 |Q| n=1 u n,k ln Pr(k|z n ; w k ) + |Q| K k=1 ûk ln(û k ), (6)
where c = stands for equality, up to an additive constant independent of optimization variables. Minimizing the last term in the second line of ( 6) has an effect opposite to the model-complexity term in our objective in (2): it reaches its minimum for perfectly balanced partitions. Therefore, our objective in ( 2) is a way to mitigate such a bias in K-means, allowing imbalanced partitions. This suggests setting λ = |Q| to compensate for the hidden class-balance term in K-means. This would make our inference hyper-parameter free. Thus, for our study, we fixed λ = |Q| in all benchmarks, without optimizing the parameter λ via validation.

Primal-Dual Block-Coordinate Descent Optimization

We derive a block coordinate descent algorithm for our formulation in [START_REF] Finn | Model-agnostic meta-learning for fast adaptation of deep networks[END_REF]. At each iteration of our algorithm, the minimization steps are closed-form, with a linear complexity in K and N . We emphasize that, in problem (2), the minimization over U could not be carried in closed-form, even for fixed prototypes W , and the simplex constraints are difficult to handle. One straightforward iterative solution to tackle (2) would be to deploy a projected gradient descent algorithm. However, as shown in the comparisons in Section 5, this strategy is computationally demanding and does not yield better classification results than the algorithm we propose below.

Primal-Dual formulation

We transform problem (2) into an equivalent optimization problem by introducing a dual variable. Let us first define the linear operator A : (R K ) N -→ R K , which maps the probabilities to the proportions of the classes, i.e. A : U → (û k ) 1≤k≤K . Then, problem (2) can be re-written as follows:

minimize U ,W 1 2 K k=1 N n=1 u n,k w k -z n 2 -λH(AU ), s.t (C)
where H is the negative entropy function on the positive orthant of R K :

H(x) = K k=1 ϕ(x k ) for x = (x k ) 1≤k≤K ∈ R K , with ϕ(t) = t ln(t) if t > 0, 0 if t = 0, +∞ otherwise. (7)
We can now appeal to the concept of Fenchel-Legendre conjugate function. For a given convex function f , its conjugate function, denoted by f * , is defined as

f * (V ) = sup U ( V , U -f (U )),
-λH(AU ) = λ inf V {H * (V ) -V , AU } . (8)
Plugging ( 8) into [START_REF] Hu | Leveraging the feature distribution in transfer-based few-shot learning[END_REF], and using the expression of H * [35, Ex. 3.21], one obtains the following minimization problem with respect to the primal variable U , the weights W and the dual variable V :

minimize U ,W ,V 1 2 K k=1 N n=1 u n,k w k -z n 2 + λ K k=1 e v k -1 -λ V , AU , s.t (C) (9)
Handling the simplex constraint To deal efficiently with the simplex constraint in (9), we add an entropic barrier on soft assignment variables u n , leading to the following modified problem:

minimize U ,W ,V 1 2 K k=1 N n=1 u n,k w k -z n 2 + λ K k=1 e v k -1 -λ V , AU + N n=1 K k=1 ϕ(u n,k ) entropic barrier , (10) s.

t (C).

The last term in [START_REF] Veilleux | Realistic evaluation of transductive few-shot learning[END_REF] acts as a barrier for imposing constraints u n ≥ 0 and, at each iteration, yields closed-form updates of both the dual variables for constraints K k=1 u n,k = 1 and assignments U . This simplifies the iterative constrained-minimization steps over U in our algorithm below through simple closed-form softmax operations (derivation details provided in supplemental material).

Block coordinate descent algorithm To minimize the cost function in [START_REF] Veilleux | Realistic evaluation of transductive few-shot learning[END_REF], we pursue a block coordinate descent approach [START_REF] Bertsekas | Nonlinear Programming, ser. Athena Scientific optimization and computation series[END_REF][START_REF] Xu | A block coordinate descent method for regularized multiconvex optimization with applications to nonnegative tensor factorization and completion[END_REF], which is guaranteed to converge; see Proposition 1 below. At each iteration, we successively minimize the objective with respect to the variables U , V , and W respectively. To do so, we use the adjoint operator of A, which we denote A * . Our algorithm is detailed in Algorithm 1.

Algorithm 1: PrimAl Dual Minimum Description LEngth (PADDLE)

Initialize W (0) as the prototypes computed on the support, and V (0) = 0. for = 1, 2, . . . , do

U ( ) = softmax - 1 2 w k -z n 2 1≤n≤N 1≤k≤K + λA * V ( -1) , v ( ) k = 1 + ln((AU ( ) ) k ), ∀k ∈ {1, . . . , K}, w ( ) k = N n=1 u ( -1) n,k z n N n=1 u ( -1)
n,k , ∀k ∈ {1, . . . , K}.

Convergence guarantees In addition to its practical advantages in terms of implementation, our algorithm benefits from the convergence guarantee described in Proposition 1.

Proposition 1. The sequence U ( ) , W ( ) , V ( ) ∈N * generated by Algorithm 1 is bounded. Moreover, any of its cluster points is a critical point to the minimization problem in [START_REF] Veilleux | Realistic evaluation of transductive few-shot learning[END_REF].

A detailed proof of Proposition 1 is provided in the supplementary material.

Experiments

Experimental details

Datasets. We deployed three datasets for few-shot classification: mini-Imagenet [START_REF] Russakovsky | ImageNet Large Scale Visual Recognition Challenge[END_REF], tiered-Imagenet [START_REF] Ren | Meta-learning for semi-supervised few-shot classification[END_REF], and i-Nat [START_REF] Wertheimer | Few-shot learning with localization in realistic settings[END_REF]. A subset of the ILSVRC-12 data [START_REF] Russakovsky | ImageNet Large Scale Visual Recognition Challenge[END_REF], mini-Imagenet is a standard few-shot benchmark, with 60, 000 color images of size 84 × 84 pixels [START_REF] Vinyals | Matching networks for one shot learning[END_REF]. It contains 100 classes, each represented with 600 images. We followed the standard split of 64 classes for base training, 16 for validation, and 20 for testing [START_REF] Ravi | Optimization as a model for few-shot learning[END_REF][START_REF] Wang | Simpleshot: Revisiting nearest-neighbor classification for few-shot learning[END_REF]. The tiered-Imagenet is another standard few-shot benchmark, which is a larger subset of ILSVRC-12, with 608 classes and a total of 779, 165 color images of size 84 × 84 pixels. We used a standard split of 351 classes for base training, 97 for validation, and 160 for testing. Finally, the more realistic and challenging dataset i-Nat has 908 classes. We follow the split from [START_REF] Wang | Simpleshot: Revisiting nearest-neighbor classification for few-shot learning[END_REF][START_REF] Ziko | Laplacian regularized few-shot learning[END_REF], with 227 ways at test-time.

Task generation

We build the few-shot tasks as follows. Let s denotes the number of shots. We constitute the support set by sampling s images according to the uniform distribution for each of the K possible classes of the test set (i.e. K = |Y test | = 20, 160, and 227 for mini-ImageNet, tiered-ImageNet and i-Nat, respectively). For the query set, we first randomly pick K eff < K classes among the K possible classes. We then randomly choose |Q| samples among the images that belong to the K eff classes but do not appear in the support set. The right-hand side of Figure 1 depicts an example of a task, with the class color coding indicating whether a support class appears in the query set or not. The results presented below are obtained by fixing |Q| = 75, as done in the literature. In the tables, we use a fixed K eff =5, but we also present the results over a larger range in Fig. 2. Given the difficulty of K-way tasks, with K large, we consider 5-, 10-and 20-shot supervision to form the support set. Following standard evaluation, we report the accuracy, averaged across 10, 000 tasks.

Hyper-parameters

We emphasize that PADDLE inference does not require any hyperparameter tuning since λ in objective ( 10) is set to the number of samples in the query set Q. This is typically not the case for other few-shot methods. For the tuning phase of those methods, we have followed the protocol of [START_REF] Veilleux | Realistic evaluation of transductive few-shot learning[END_REF], which for each dataset, uses a 5-way, 5-shot scenario on the corresponding validation set. For i-Nat, because no validation split is provided, we reuse the hyper-parameters obtained from tiered-ImageNet.

Feature extraction To ensure the fairest comparison of methods, we use the pretrained checkpoints provided by the authors of [START_REF] Boudiaf | Transductive information maximization for few-shot learning[END_REF]. The models are trained on D base via a standard cross-entropy minimization with label smoothing. The label smoothing parameter is set to 0.1, for 90 epochs, using a learning rate initialized to 0. Results are averaged across 10,000 tasks. Results marked with '-' were intractable to obtain.

Results

Main results

We compare the performances of our method with state-of-the-art few-shot methods. Our first experimental setting consists in fixing the number of effective classes K eff to 5. In Table 1, we evaluate the accuracy on mini and tiered-ImageNet for 5, 10, and 20 shots, while results on i-Nat are displayed in Table 2. Note that the i-Nat dataset comes with a unique support set for all tasks, with a varying number of shots (labeled samples) per class. Therefore, for i-Nat, we present the results separately. In the second experiment displayed in Figure 2 and 3, we plot the accuracy as a function of K eff on mini, tiered-ImageNet, and i-Nat for 5, 10, and 20 shots. Both experiments show that PADDLE is very competitive with the state-of-the-art methods on the proposed new practical few-shot setting. More importantly, the gap between PADDLE and the second best algorithm increases significantly with (i) the number of shots and (ii) the number of possible classes K, when going from the dataset mini, which has a small number of test classes (K = 20), to tiered (K = 160) and finally to i-Nat (K = 227). From the first point, one may conclude that the proposed MDL formulation benefits better from additional support supervision. We hypothesize to this is due to its versatility, as it does not encode strong assumptions on the label statistics of the query set. Also, the second point could be explained along the same line. For instance, methods relying on a class-balance assumptions, such as TIM-GD or PT-MAP, observe a drastic drop in accuracy as K increases, even falling below the inductive Baseline. This is expected because large values of K corresponds to highly imbalanced classification problems for the query sets. For example, K eff =5 effective classes on tiered corresponds to a relatively more drastic form of class imbalance than on mini, as it would mean only 5/160 classes are represented in Q, versus 5/20. In contrast, PADDLE can deal with strongly imbalanced situations. Note that BD-CSPN [START_REF] Liu | Prototype rectification for few-shot learning[END_REF] is the only baseline whose performance increases with K eff . This method is also among the most affected when the class overlap between the query and support sets decreases, i.e., when K eff gets smaller (left side of the plots in Figure 2) and/or K gets bigger (e.g. tiered-ImageNet, with K = 160). This behaviour might be due to the fact that BD-CSPN encodes a strong prior, assuming the support and query classes match perfectly. Baseline [START_REF] Chen | A closer look at few-shot classification[END_REF] 58.2 BD-CSPN [START_REF] Liu | Prototype rectification for few-shot learning[END_REF] 57.6 PT-MAP [START_REF] Hu | Leveraging the feature distribution in transfer-based few-shot learning[END_REF] 6.8 LaplacianShot [START_REF] Ziko | Laplacian regularized few-shot learning[END_REF] 43.3 TIM [START_REF] Boudiaf | Transductive information maximization for few-shot learning[END_REF] 37.5 α-TIM [START_REF] Veilleux | Realistic evaluation of transductive few-shot learning[END_REF] 66.7 PADDLE (ours) 84.3 Ablation on the objective We hereby ablate on the importance of the partition-complexity term in Eq. [START_REF] Finn | Model-agnostic meta-learning for fast adaptation of deep networks[END_REF]. Removing this high-order term yields a partially-supervised version of K-means, which can be optimized effortlessly through iterative closed-form alternating steps. We provide the comparison between this approach (without the partition-complexity term) and PADDLE (with the partition complexity term) in Table 3. On can observe that PADDLE systematically outperforms its K-Means counterpart, with absolute differences in accuracy reaching up to 30% + on challenging scenarios. Table 3: Importance of the partition-complexity term. Results are computed using a ResNet-18 and K eff =5, and averaged across 10,000 tasks. Without the partition-complexity term in Eq. ( 2), the algorithm has little incentive to be parsimonious in its choice of the effective classes. As a matter of fact, removing this term reduces Problem (2) to a partially-supervised K-means algorithm, notoriously known to encourage balanced solutions across the classes.

Ablation on the optimization procedure As a second ablation, we propose to compare our proposed Alternating Minimization Algorithm 2 to a straightforward first-order approach. More precisely, as a comparison, we directly optimize Eq. ( 2) through Projected Gradient Descent, denoted as PGD , where a simplex projection step [START_REF] Bauckhage | NumPy/SciPy recipes for data science: Projections onto the standard simplex[END_REF] takes place after every iteration to ensure that the simplex constraints on U remains satisfied. For PGD , we use Adam [START_REF] Kingma | Adam: A method for stochastic optimization[END_REF] with α = 0.001. We found this to be the highest learning rate that leads to convergence. We also add the convergence time of the second best competing method, α-TIM. Results are provided in Fig. 4. PGD is much slower than PADDLE , with a ratio between run times neighboring an order of magnitude. Additionally, as hinted from the last points of the PGD curve, oscillations seem to occur, indicating that a more sophisticated learning rate policy would be necessary to achieve a clear convergence.

Conclusion and limitations

We presented a practical few-shot setting where the number of candidate classes can be much larger than the number of classes that appear effectively in the query set. Our setting is an instance of highly imbalanced classification, with large numbers of ways and potentially irrelevant supervision from the support set. We observed much higher gaps across transductive few-shot methods, some of which fall below the simple inductive baseline. As a solution, we introduced PADDLE , which casts this challenge as a partially-supervised MDL partitioning problem, interpreting the number of unique classes found as a measure of model complexity. PADDLE is hyperparameter-free, and remains competitive over various state-of-the-art methods, settings and datasets, without the need for any tuning. However, we do not advocate our method as the one-fits-all, ultimate solution to the extremely challenging few-shot problem. For instance, PADDLE would not be the best performing method in situations where K eff would approach K. In applications where one has knowledge about the label statistics of the query set (e.g. class balance), other methods encoding such a knowledge could be more appropriate. More generally, our setting can be seen as a particular case of class-distribution shift between the support and query sets. Interesting future works could complement our current study by overlaying other forms of shifts, e.g. feature shifts [START_REF] Guo | A broader study of cross-domain few-shot learning[END_REF][START_REF] Bennequin | Bridging few-shot learning and adaptation: New challenges of support-query shift[END_REF], to our setting.

Supplementary material A Closed-form updates of the assignment variables

In this section, we provide more details on the derivation of the closed-form update of variable U at each iteration. Let F be the defined as the cost function in [START_REF] Veilleux | Realistic evaluation of transductive few-shot learning[END_REF] and let ∂F un (U , W , V ) denote the Moreau subdifferential of F at (U , W , V ) with respect to variable u n . We define ψ as

(∀x = (x k ) 1≤k≤K ∈ R K ) ψ(x) =      K k=1 x k ln(x k ) - x 2 k 2 if x ∈ ∆ K , +∞ otherwise. (11)
It is well known that the proximity operator of ψ (see [START_REF] Bauschke | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF]Chap. 24] for a definition) is the softmax operator [START_REF] Combettes | Deep neural network structures solving variational inequalities[END_REF]Ex. 2.23].

At each step of the algorithm, u n is updated according to:

0 ∈ ∂F un (U , W , V ) ⇐⇒ 0 ∈ 1 2 w k -z n 2 1≤k≤K -λ[A * V ] n + u n + ∂ ψ (u n ), ⇐⇒ - 1 2 w k -z n 2 1≤k≤K + λ[A * V ] n -u n ∈ ∂ ψ (u n ), ⇐⇒ u n = softmax - 1 2 w k -z n 2 1≤k≤K + λ[A * V ] n , (12)
where we used the definition of the proximity operator [34, Eq. 24.2] to obtain [START_REF] Boudiaf | Few-shot segmentation without meta-learning: A good transductive inference is all you need?[END_REF]. We thus retrieve the update in Algorithm 1.

B Proof of Proposition 1

Our proof relies on the convergence result established in [START_REF] Tseng | Convergence of a block coordinate descent method for nondifferentiable minimization[END_REF]. Given a convex set X, we denote ι X the indicator function of X, i.e. ι X (x) = 0 if x ∈ X, ι X (x) = +∞ otherwise. We rewrite problem 10 as the minimization of the following cost:

F (U , W , V ) = 1 2 K k=1 N n=1 u n,k w k -z n 2 + λ K k=1 e v k -1 -λ V , (AU + 1 K ) + N n=1 K k=1 ϕ(u n,k ) + ι C (U ), (13)
where we have introduced an additional parameter > 0, the role of which will become clearer in the rest of the proof. The optimum of the cost function

F (U , W , •) for given U ∈ C and W ∈ (R d ) K is reached when V = 1 K + ln(AU + 1 K ) ∈ V = [1 + ln , 1 + ln(1 + )] K . (14)
Thus, minimizing F is actually equivalent to minimizing

F (U , W , V ) = 1 2 K k=1 N n=1 u n,k w k -z n 2 + λ K k=1 e v k -1 -λ V , (AU + 1 K ) + N n=1 K k=1 ϕ(u n,k ) + ι C (U ) + ι V (V ). (15)
The following algorithm for minimizing F turns out to be a simple modified version of PADDLE (see Algorithm 1):

The latter function is strictly convex, lower-semicontinuous, and coercive, which concludes the proof. Note that, since

v k → λ e v k -1 -v k ([AU ( +1) ] k + ) (19)
is decreasing on ] -∞, 1 + ln([AU ( +1) ] k + )] and increasing on [1 + ln([AU ( +1) ] k + ), +∞[, the resulting cluster points are also coordinatewise minimizers of F . In summary, PADDLE can be understood as the limit case of Algorithm 2 when goes to zero. This simplification is justified by the fact that can be chosen arbitrarily small and that we did not observe any change in practical behaviour of the proposed algorithm by setting = 0.

C Label cost relaxation

The plot in Figure 5 illustrates in the case K = 2 how our model-complexity term in (2) could be viewed as a continuous relaxation of the discrete label cost function defined in (3).

D Plots obtained using WRN backbone

In Figure 6, we provide additional comparisons of PADDLE with state-of-the-art methods using a WRN28-10 network. We report the accuracy obtained for each method as a function of K eff . These plots point to the same conclusions drawn in Section 5.

E About the hyper-parameter in our method

As discussed in Section 3, PADDLE does not require parameter tuning. In Figure 7, we investigate the optimal value of parameter λ in (10) as a function of the size of the query set, for 3 different values of K eff . We observe that the optimal value of λ increases linearly with |Q|. As it could be expected, the higher the level of class imbalance (K eff = 2), the higher the optimal value of λ (w.r.t. its theoretical value). On the contrary, when the query is better balanced (K eff = 10), the optimal value of λ is slightly under its theoretical value. However, Figure 8 shows that the gap of performance when using the theoretical value of λ instead of the optimal one, is only of the order of a few percents.

Figure 1 :

 1 Figure 1: Overview of the proposed framework: Base training, PADDLE inference and task generation.The example depicted in the right-hand side illustrates how the support and query classes do not match exactly, unlike in standard few-shot settings. The support set includes "distraction" classes that may not actually be present in the query set, e.g. classes "Alligator", "Peacock" and "Butterfly". All possible classes (K classes) are represented in the support set, but only an unknown subset among these K possible classes (K eff effective classes) appear in the query set, with K eff K.

Figure 2 :

 2 Figure 2: Evolution of the accuracy as a function of K eff . Each row represents a dataset, and each column a fixed number of shots. All methods use the same ResNet-18 network. Results are averaged across 10,000 tasks.

Figure 3 :

 3 Figure 3: Similarly to Fig. 2, we plot the performances of the methods as functions of K eff on i-Nat with 227-ways tasks.

Figure 4 :

 4 Figure 4: Comparison of PADDLE and PGD in terms of convergence speed (Left) The chosen criterion as a function of the elapsed time on a randomly chosen 20-shot, K eff =5 task sampled from tiered-ImageNet. Both methods are run on the same machine, with markers displayed every 100 iterations. (Right) time to convergence and accuracy, after reaching a value of 10 -6 on the criterion.

Figure 5 :

 5 Figure 5: Label cost as a function of û1 and our proposed relaxation û1 → -û 1 ln(û 1 ) -(1 -û1 ) ln(1 -û1 ).

Figure 6 :PADDLEFigure 7 :PADDLEFigure 8 :

 678 Figure 6: Evolution of the accuracy as a function of K eff . Each row represents a dataset, and each column a fixed number of shots. All methods use the same WRN28-10 network. Results are averaged across 10,000 tasks.

  Evaluation Evaluation is carried out over few-shot tasks, each containing samples from D test = {x n , y n } |Dtest| n=1 , where y n ∈ {0, 1} |Ytest| , with constraint Y base ∩ Y test = ∅, i.e., the test and base classes are distinct. Each task includes a labelled support set S = {x n , y n } n∈I S and an unlabelled query set Q = {x n } n∈I Q , both containing examples from classes in Y test . Using the feature extractor f φ trained on base data, the goal is to predict the classes of the unlabeled samples in Q for each few-shot task, independently of the other tasks.

  whereV is a dual variable[START_REF] Bauschke | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF] Def. 13.1] and •, • denotes the Euclidean scalar product on R K×N . Since function H is a proper, lower semicontinuous, convex function, according to the Fenchel-Moreau theorem[START_REF] Bauschke | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF] Thm. 13.37], the equality H = (H

* ) * holds. Hence, it follows that

Table 1 :

 1 1 and divided by 10 at epochs 45 and 66. We use batch sizes of 256 for ResNet-18 and of 128 for WRN28-10. The images are resized to 84 × 84 pixels, both at training and evaluation time. Color jittering, random croping, and random horizontal flipping augmentations are applied during training. Comparisons of state-of-the-art methods on mini-Imagenet and tiered-Imagenet, using the tasks generation process described in Sec. 2 with K eff =5. The metric is accuracy (in percentage).

	Method	Backbone	mini-ImageNet (K =20)	tiered-ImageNet (K =160)
			5-shot 10-shot 20-shot 5-shot 10-shot 20-shot
	Baseline [24]		54.1	60.7	65.8	29.1	35.7	39.5
	LR+ICI [15]		54.7	62.0	67.2	-	-	-
	BD-CSPN [8]		49.6	54.1	55.6	24.0	26.7	26.0
	PT-MAP [7]	ResNet-18	25.8	27.4	29.0	4.3	5.1	5.9
	LaplacianShot [4]		60.4	65.2	68.4	34.7	36.8	39.1
	TIM [11]		66.4	69.4	70.9	26.6	26.7	25.5
	α-TIM [10]		63.5	67.4	71.7	38.7	44.2	48.4
	PADDLE (ours)		63.2	73.3	80.0	45.8	62.5	72.0
	Baseline [24]		58.0	64.8	69.5	31.8	37.0	42.1
	LR+ICI [15]		57.2	64.3	70.9	-	-	-
	BD-CSPN [8] PT-MAP [7]	WRN28-10	51.2 26.3	55.5 27.9	58.4 29.4	23.7 4.4	24.9 5.0	23.8 5.7
	LaplacianShot [4]		64.9	65.3	70.9	28.1	38.0	45.2
	TIM [11]		71.9	75.2	76.1	34.1	34.1	34.5
	α-TIM [10]		68.3	72.5	75.8	41.9	46.1	51.8
	PADDLE (ours)		62.5	72.8	79.3	46.4	60.3	71.1

Table 2 :

 2 Similarly to Table 1, results are provided with a fixed K eff =5 on i-Nat with 227-ways tasks.

Note that the transductive setting is different from semi-supervised few-shot learning[START_REF] Ren | Meta-learning for semi-supervised few-shot classification[END_REF] that uses extra data. The only difference between inductive and tranductive inference is that predictions are made jointly for the query set of a task, rather than one sample at a time.

In the general context of clustering, the name probabilistic K-means was first coined by[START_REF] Kearns | An information-theoretic analysis of hard and soft assignment methods for clustering[END_REF].

* E. Chouzenoux acknowledges support from the European Research Council Starting Grant MAJORIS ERC-2019-STG-850925. † The work of J.-C. Pesquet is supported by the ANR Chair in AI BRIDGEABLE. ‡ The work of I. Ben Ayed is supported by the DATAIA Institute, and is part of his sabbatical-leave visit at the Université Paris-Saclay.

Our code is publicly available at https://github.com/SegoleneMartin/PADDLE.

Algorithm 2: Alternating algorithm for minimizing F Initialize W (0) as the prototypes computed on the support, and V (0) = 0. for = 1, 2, . . . , do

n,k , ∀k ∈ {1, . . . , K}.

According to [START_REF] Tseng | Convergence of a block coordinate descent method for nondifferentiable minimization[END_REF]Thm 4.1], if the following assumptions are satisfied:

3. At each iteration , the partial functions F (•, W ( ) , V ( ) ), F (U ( +1) , •, V ( ) ) and F (U ( +1) , W ( +1) , •) admit a unique minimizer, then the sequence generated by the algorithm is bounded and every of its cluster points is a coordinatewise minimizer of F . We now show that the above assumptions hold.

1. Let us show that F is coercive. We derive a lower bound on F using the Cauchy-Schwarz inequality:

Since the functions U → AU and U → N n=1 K k=1 ϕ(u n,k ) are continuous on the compact set C, there exist constants µ and θ such that

The lower bound obtained in [START_REF] Shen | Re-ranking for image retrieval and transductive few-shot classification[END_REF] is separable in (U , W , V ). The term with respect to variable W is coercive when, for every k ∈ {1, . . . , K}, there exists n ∈ {|Q| + 1, . . . , N } such that y n,k > 0. In other words, it is coercive if the support set includes at least one example of each class, which is a reasonable assumption. The terms with respect to variables U and V are clearly coercive too. Hence, the cost function F is coercive. Finally, since F is lower semi-continuous, condition 1. is satisfied. 2. The continuity of F on C × R k×d × V is clear. 3. Let ∈ N * . We already proved in Appendix A that the partial function with respect to variable U has a unique minimizer. It follows from the same arguments as above that the partial function with respect to W is strictly convex, continuous, and coercive as soon as the support set contains at least one example of each class. Hence, it admits a unique minimizer.

Regarding the partial function with respect to variable V , we first remark that given the definition of the softmax operator, AU ( +1) is necessarily strictly positive component-wise.

Up to some additive term independent of V , the partial function reads


Céline Béchaud

Stéphane Mélen

Didier Lasseux
email: didier.lasseux@u-bordeaux.fr

Michel Quintard

Charles-Henri H Bruneau

C Bã

S Mã Elen C

Stability analysis of a pressure swing adsorption process

Keywords: Adsorption, Dispersion, Porous media, Stability, Transport processes, Zeolites

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Pressure swing adsorption processes-referred to as PSA in this work-are widely used in a variety of industrial applications like air drying, gas puriÿcation, solvent recovery among the principals. For oxygen production from air in particular, this technique o ers interesting advantages in terms of energy consumption, as compared to cryogenic processes for instance, and allows on-site installation due to a relatively light infrastructure. Basically, the technique consists in adsorbing nitrogen of air on reactive sites within a porous bed made of zeolite to produce oxygen at the bed outlet.

Originally, PSA is based on the so-called [START_REF] Skarstrom | Use of adsorption phenomena in automatic plant-type gas analysers[END_REF] cycle, even though many other cycles were derived afterwards for performance purposes [START_REF] Turnock | Separation of nitrogen and methane via periodic adsorption[END_REF][START_REF] Kowler | The optimal control of a periodic adsorber[END_REF]. The Skarstrom cycle involves two main steps: an adsorption step during which oxygen is actually produced, and a depressurization one which corresponds to a regeneration of the bed.

The overall process is characterized by its performance which is usually estimated with the aid of two main parameters: product recovery and adsorbent productivity. The former represents the concentration ratio between the inlet and outlet of the non-adsorbed componentoxygen for instance-while the latter is the amount of the non-adsorbed component produced per unit of adsorbent. Many studies on the PSA process have been carried out in order to improve its performance. For instance, [START_REF] Turnock | Separation of nitrogen and methane via periodic adsorption[END_REF] and [START_REF] Kowler | The optimal control of a periodic adsorber[END_REF] developed a variation of the Skarstrom one-column process (called pressure swing parametric pumping process). This was designed to separate a mixture of nitrogen and methane. They focused their attention on the feed pressure wave-form, and observed experimentally that the optimum feed pressure wave-form was a square one. Considerable improvements of this parametric pumping process, in the case of oxygen and nitrogen separation from air, was obtained later by [START_REF] Jones | Pressure swing parametric pumping-a new adsorption process[END_REF] by modifying the time steps in the cycle. Another idea to increase PSA performance has been proposed. It consists in coupling the process with other gas separation techniques. For instance, Vaporciyan andKadlec (1987, 1989), [START_REF] Lu | Pressure swing adsorption reactors: Simulation of three-step one-bed process[END_REF] combined chemical reaction to PSA, and [START_REF] Bhaumik | Hollow-ÿber membrane-based rapid pressure swing adsorption[END_REF] integrated the best features of membrane contacting, gas liquid absorption and PSA to give rise to the so-called rapid pressure swing absorption (RAPSAB). Other investigations were performed on the PSA performance. E ciency of various adsorbents, like carbon molecular sieves (CMS) and zeolite, was compared by [START_REF] Farooq | A comparison of linear driving force and pore di usion models for a pressure swing adsorption bulk separation process[END_REF]. They found that CMS lead to a higher product recovery whereas zeolite productivity is better. For this reason, most of the PSA processes are designed with zeolite as adsorbent material. Later, [START_REF] Lemco | E ect of regeneration pressure level in kinetically controlled pressure swing adsorption[END_REF] and [START_REF] Shirley | Novel pressurization methods in pressure swing adsorption systems for the generation of high purity gas[END_REF] stressed the importance of the regeneration pressure on the performance. Their conclusion was that a better performance is achieved when pressurization is done by feed gas and product backÿll. E ects of parameters like "dead" zones, adsorption isotherms, incomplete pressurization and blow-down, as well as feed composition, isothermal or adiabatic operations on the separation performance were investigated by [START_REF] Lu | Simulation of a three-step one-column pressure swing adsorption process[END_REF]. They showed that the number of cycles required to achieve the cyclic steady state is shorter in the isothermal case than in the adiabatic one. A recent study of [START_REF] Shirley | High-purity nitrogen by pressure-swing adsorption[END_REF] highlighted the e ect of cycle time on performance as well as on purity of the product. It was found that performance improves with increasing cycle time especially for high purity whereas purity is independent of cycle time for low purity regimes.

In existing works, two main approaches have been followed to model the process. On the one hand, the pore di usion model was introduced to represent the intra-particle di usional mass transfer resistance. This model was used by [START_REF] Yang | Gas separation by pressure swing adsorption: A pore-di usion model for bulk separation[END_REF], [START_REF] Sun | An improved ÿnite di erence method for ÿxed-bed multicomponent sorption[END_REF] and [START_REF] Lu | Simulation of a three-step one-column pressure swing adsorption process[END_REF] to predict concentration evolution during the di erent steps. On the other hand, a more simpliÿed model, referred to as the linear driving force (LDF) model, has been used. It is based on an approximated description of the ow within the pores, and does not take explicitly into account molecular di usion [START_REF] Carter | The pressure swing adsorption drying of compressed air[END_REF]. Comparisons between these two models and experimental data were performed by [START_REF] Farooq | A comparison of linear driving force and pore di usion models for a pressure swing adsorption bulk separation process[END_REF]. Their results indicated that the pore di usion model provides a better agreement with actual data than the LDF model, in spite of a significantly longer computational time. Di erences observed between model predictions and experimental data may result from simpliÿcations involved in the physical model. They may also be explained by additional mechanisms that are not taken into account in these models. Among the most important ones, heterogeneity e ects or development of hydrodynamic instabilities may a ect dramatically the ow. In this paper, we consider the possible occurrence of instabilities, and we investigate the in uence of various physical parameters involved in the process on the development of such instabilities.

Hydrodynamic instabilities have been reported for a long time in the literature in the case of passive dispersion in porous media (see a review by [START_REF] Homsy | Viscous ÿngering in porous media[END_REF], or in the case of the development of a thermal di usion zone (see [START_REF] Quintard | Stabilità e d'une zone instationnaire de di usion thermique dans un milieu poreux: Analyse linà eaire, modà elisation numà erique[END_REF]. Interaction of these instability mechanisms with heterogeneities in the porous bed has also been investigated [START_REF] Homsy | Viscous ÿngering in porous media[END_REF][START_REF] Tan | Viscous ÿngering with permeability heterogeneity[END_REF]. However, occurrence of instabilities during dispersion and adsorption, and for the case of compressible gases, has not been studied so far. In the present work, a stability analysis was performed on the basis of a pore di usion model used to describe the transport of the two chemical species-oxygen and nitrogen. The regular ow corresponds to the 1D solution of the governing equations. Since our objective is to develop the analysis in the case of PSA processes, the resulting 1D base ow is time-dependent, periodic in time, with a complex evolution due to the periodic change in the boundary conditions. To investigate the stability of this time-dependent solution, a linear amplitude equation associated with 2D inlet perturbations was developed [START_REF] Quintard | Stabilità e d'une zone instationnaire de di usion thermique dans un milieu poreux: Analyse linà eaire, modà elisation numà erique[END_REF]. A numerical model was derived to solve for this amplitude equation. Time evolution of the disturbance amplitude indicates whether instabilities can develop or not. Results allow us to discuss the e ect of the perturbation wavelength and amplitude on the development of instabilities.

Physical problem

The system under consideration was designed to produce oxygen from air (mixture of 78% nitrogen and 21% oxygen). Air ows within three columns containing a bed of porous zeolite beads. Zeolite adsorb nitrogen more than oxygen. The air-ow cycle under study derives from the Skarstrom cycle, and corresponds to three main steps that are detailed below. To ensure a continuous production [START_REF] Yang | Gas separation by pressure swing adsorption: A pore-di usion model for bulk separation[END_REF], all three columns are used simultaneously, each of them undergoing one cycle step (see Fig. 1).

The three main steps of a cycle correspond to:

(1) adsorption during two distinct periods: (I) a feeding period: air is injected and adsorption of nitrogen results in oxygen production; (II) a delay period during which adsorption continues while depressurization begins; (2) purge:

(III) during this period, the bed is depressurized with the help of a vacuum pump; (3) pressurization during which one can observe:

(IV) an elution period which corresponds to a pressurization with gas coming from the bed running in period (II). A consequence of this design is that elution period has the same time length as period II; (V) a compression period to regenerate the bed using oxygen produced in period (I).

Production of oxygen occurs during the two ÿrst periods, while the three other ones are designed to regenerate the bed in order to restore the adsorption capacity of zeolite. The overall process is represented schematically in Fig. 2 where we have represented the evolution of the entrance pressure. It is important to notice that all the three columns are coupled so that gas production at the outlet of one column is used as inlet for another one.

Physical mechanisms in this problem can be considered at many di erent scales: the pore scale which corresponds to pores within the porous zeolite beads; the macro-pore scale, with a characteristic length of beads diameter; the scale of few beads corresponding to an homogeneous pack which is often referred to as Darcy-scale in the porous media literature; and the mega-scale corresponding to the entire column with packing heterogeneities. In the following, the solution to the stability problem is proposed at the macroscopic-Darcy-scale and equations at this scale are taken from the literature. The problem of deriving these macroscopic equations from some averaging technique applied to the pore-scale equations is beyond the scope of this paper (see, for instance, [START_REF] Whitaker | Flow in porous media. I: A theoretical derivation of Darcy's law[END_REF][START_REF] Hassanizadeh | General conservation equations for multi-phase systems: 2. Mass, momentum, energy and entropy equations[END_REF][START_REF] Quintard | Convection, dispersion and interfacial transport of contaminants: Homogeneous porous media[END_REF], and they are taken without further discussion. For the initial boundary-value problem, we assume that physical parameters of the porous bed, like porosity, dispersion coe cients, heat capacities and heat conduction coe cients are constant. Since the column diameter is much smaller than its length, the ow problem is considered as 1D and is solved on the homogeneous system schematically depicted in Fig. 3. Because gas phase density is small, gravity is neglected since this type of force is assumed to play a minor role compared to adsorption forces.

The gas phase, or ÿ-phase, is a binary mixture of oxygen, ÿ 1 , and nitrogen, ÿ 2 , while the solid phase,

, is assumed to be rigid and immobile. The initial boundary-value problem is given by mass, momentum and energy balance equations, associated with boundary conditions corresponding to the cycle under study. The set of equations are given below.

Mass balance

Mass balance for component i (i = ÿ 1 or ÿ 2 ) in the gas mixture is a convection=dispersion equation with additional terms taking into account gas compressibility and adsorption. For the oxygen component, this writes

@ ÿ @t C ÿ1 + ÿ @C ÿ1 @t compression + (1 -) @f(C ÿ1 ) @t adsorption + • ( ÿ V ÿ C ÿ1 ) convection = • ( ÿ D * ÿ • C ÿ1 ) dispersion ; (1)
where C ÿ1 is the macroscopic oxygen mass fraction in the gas phase deÿned by

C ÿi = m i m ; (2)
m i being the mass of component i (i = 1; 2) in a mass m of the mixture; V ÿ and D * ÿ are the macroscopic velocity and dispersion tensor; is the porosity of the medium and f represents the adsorption isotherm deÿned below. Rather than writing the equivalent equation for the nitrogen component, it is more convenient to write the total mass balance that results from the sum of Eq. ( 1) for both components to obtain

@ ÿ @t compressibility + • ( ÿ V ÿ ) convection = -(1 -) @ C 1 @t + @ C 2 @t adsorption ; (3)
where C 1 and C 2 denote, respectively, the mass of oxygen and nitrogen adsorbed per unit mass of the porous medium.

Adsorption is a physico-chemical phenomenon that depends in a nonlinear way on pressure, temperature and species concentration in the uid phase. For a given adsorbent, this relationship is generally given under the form of a so-called isotherm, and in the present work, a Langmuir relationship was used. This kind of relationship has proved to be very accurate for the air=zeolite system at low pressure [START_REF] Van Tassel | New lattice model for adsorption of small molecules in zeolite micropores[END_REF]. In addition, we assumed that the macroscopic concentration is equal to the pore-scale equilibrium concentration. This assumption corresponds to a local chemical equilibrium assumption and was discussed in some details by Quintard andWhitaker (1998), and[START_REF] Ahmadi | Transport in chemically and mechanically heterogeneous porous media V: Two-equation model for solute transport with adsorption[END_REF]Whitaker (1998). As a consequence the following Langmuir isotherm was used, which, for component i reads

C i = f(C ÿi ) = i exp( h i =RT )C m ÿi P ÿ 1 + Ä i exp( h i =RT )C m ÿi P ÿ : (4)
In this isotherm, C m ÿi represents the molecular fraction of component i and is given by

C m ÿi = n i n m ; ( 5
)
where n i is the number of molecules of component i contained in a mixture of n m molecules; h i represents the molar enthalpy of adsorption; i and Ä i are system dependent constants and R is the ideal gas constant.

Momentum balance

Provided the pore Reynolds number is less than 10, it is usually admitted that a valid model for momentum conservation in a homogeneous porous medium is the classical Darcy's law. In our case, the Reynolds number estimated from the mean pore diameter and interstitial velocity is below this critical value and consequently, since gravity is neglected, the following relationship was used for the ÿ-phase momentum balance:

V ÿ = - 1 ÿ K • P ÿ ; ( 6
)
where K is the permeability tensor deÿned on a homogeneous representative sample of the bed.

Energy balance

Following some ideas put forth for the local chemical equilibrium assumption, local thermal equilibrium between the gas phase and the solid phase was further assumed [START_REF] Quintard | Two-medium treatment of heat transfer in porous media: Numerical results for e ective properties[END_REF]. This yields a unique macroscopic equation to describe the temperature evolution within the uid-solid system which reads ( C p ) * @T @t accumulation term

+ • (( C p ) ÿ V ÿ T ) convection = • (K * • T ) heat di usion + (1 -) @ C 1 @t h 1 + @ C 2 @t h 2 adsorption + @P ÿ @t + V ÿ • P ÿ pressure reduction : (7)
In the above equation, K * and ( C p ) * represent the equivalent heat conductivity tensor and the equivalent volumetric heat capacity of the uid=solid system, respectively. This latter property is deÿned by a mixing law of the following type:

* = ÿ + (1 -) : (8)
To close the problem, one needs a constitutive equation for the gas phase. For simplicity, an ideal gas law was used to deÿne gas density. According to the deÿnitions of the mass fractions of the two components, one obtains the density of the mixture as

ÿ = P ÿ [C ÿ1 =M ÿ1 + (1 -C ÿ1 )=M ÿ2 ]RT ; (9)
where M ÿi is the molar mass of component i.

Initial and boundary conditions

Boundary conditions at entrance, e , and exit, s , of one column are those corresponding to the successive periods over a complete cycle. For sake of clarity, they are listed below for each period:

Period I (feeding):

(B:C:1) T (t) = T e on e ; (10)
(B:C:2) C ÿ1 (t) = C ÿ1 e on e ; (11)
(B:C:3) P ÿ (t) = P ÿe on e ;

(12) (B:C:4) P ÿ (t) = P ÿs on s ;

(13)

(B:C:5) n • K * • T (t) = 0 on s ; (14) (B:C:6) n • D * ÿ • C ÿ1 (t) = 0 on s : (15)
Period II (delay):

(B:C:7) V ÿ (t) = 0 on e ; (16)
(B:C:8) n • K * • T (t) = 0 on e ; (17)
(B:C:9) n • D * ÿ • C ÿ1 (t) = 0 on e ; (18)
(B:C:10) P ÿ (t) = P ÿs (t) on s ;

(B:C:11) n • K * • T (t) = 0 on s ; (20) (B:C:12) n • D * ÿ • C ÿ1 (t) = 0 on s : (19)
During this period, P ÿs (t) corresponds to a ÿxed linear pressure evolution. Values of V ÿ ; C ÿ1 and T over time at the outlet are later used in period IV.

Period III (depressurization):

(B:C:13) P ÿ (t) = P ÿe (t) on e ;

(B:C:14) n • K * • T (t) = 0 on e ; ( (22)
) (B:C:15) n • D * ÿ • C ÿ1 (t) = 0 on e ; (24) (B:C:16) V ÿ (t) = 0 on s ; (25) (B:C:17) n • K * • T (t) = 0 on s ; (26) (B:C:18) n • D * ÿ • C ÿ1 (t) = 0 on s : 23
During this depressurization period, P ÿe (t) is an imposed linear decreasing function of time.

Period IV (elution):

(B:C:19) P ÿ (t) = P ÿe (t) on e ;

(28)

(B:C:20) n • K * • T (t) = 0 on e ; (29)
(B:C:21) n • D * ÿ • C ÿ1 (t) = 0 on e ; (30)
(B:C:22) V ÿ (t) = 0 on s ; (31) (B:C:23) T (t) = T s (t) on s ; (32)
(B:C:24) C ÿ1 (t) = C ÿ1 s (t) on s : (33)
During this period, P ÿe (t) is an imposed linear increasing function of time and conditions at any time on s are those corresponding to period II.

Period V (pressurization):

(B:C:25) V ÿ (t) = 0 on e ; (34)
(B:C:26) n • K * • T (t) = 0 on e ; (35)
(B:C:27) n • D * ÿ • C ÿ1 (t) = 0 on e ; (36)
(B:C:28) P ÿ (t) = P ÿs (t) on s ;

(37)

(B:C:29) T (t) = T initial on s ; (38)
(B:C:30) C ÿ1 (t) = C ÿ1 initial on s : (39)
During this ÿnal period, P ÿs (t) is a ÿxed linear pressure evolution, whereas oxygen concentration as well as temperature on s are kept constant and equal to their initial values at the beginning of the cycle.

This initial boundary-value problem was solved numerically and a stability analysis with respect to pressure or concentration perturbations was performed. Methodology and numerical procedures used to compute the regular solution on the ÿ phase pressure, P 0 ÿ , temperature, T , oxygen concentration, C 0 ÿ1 , and velocity, V 0 ÿ , as well as the corresponding disturbances are now presented.

Hydrodynamic stability

Our study of the hydrodynamic stability of the PSA process under consideration is inspired by techniques developed for the case of miscible displacements [START_REF] Wooding | The stability of an interface between miscible uids in a porous medium[END_REF][START_REF] Bachu | Stability of displacement of a cold uid by a hot uid in a porous medium[END_REF][START_REF] Bertin | Stabilità e d'un à ecoulement miscible radial en milieu poreux: à Etude thà eorique et expà erimentale[END_REF][START_REF] Tan | Stability of miscible displacements in porous media: Rectilinear ow[END_REF][START_REF] Homsy | Viscous ÿngering in porous media[END_REF][START_REF] Yortsos | Dispersion driven instability in miscible displacement in porous media[END_REF][START_REF] Prouvost | Stability criteria for the design of graded polymer bu ers[END_REF] or transient heat diffusion [START_REF] Gresho | The stability of a uid layer subjected to a step change in temperature: Transient vs. frozen time analyses[END_REF][START_REF] Homsy | Global stability of time-dependent ows: Impulsively heated or cooled uid layers[END_REF][START_REF] Quintard | Instabilità es de zones de di usion thermique instationnaires en milieu poreux[END_REF].

A very popular technique is the classical bifurcation theory-or linear stability analysis- [START_REF] Iooss | Elementary stability and bifurcation theory[END_REF] which is well adapted to the case of autonomous regular solution. In fact, if W 0 denotes the regular solution vector or base state, and w the associated disturbance, then the equations for the disturbance can be linearized in the neighborhood of the regular solution and formally written as @w(x; t) @t = L(w(x; t); W 0 (x));

(40

)
where L is the linear operator and W 0 the regular solution which is a function of the spatial coordinates, x, only.

The stability analysis can be performed by studying the behavior of the operator L with respect to an elementary perturbation of the form:

w(x; t) = e t '(x); (41)
which implies

'(x) = L('(x); W 0 ); ( 42
)
where and ' are, respectively, eigenvalues and eigenvectors of the operator L. Marginal stability is obtained for values of the physical parameters such that = 0. In our case, however, the regular solution (we also used the term base ow) is time dependent in a complex non-linear fashion and the resulting disturbance equations lead to a non-autonomous system. A non-autonomous system is characterized by the fact that coe cients in the disturbance equations are time dependent while the regular solution is a non-stationary one. For this kind of system, more complex approaches are required [START_REF] Homsy | Global stability of time-dependent ows: Impulsively heated or cooled uid layers[END_REF][START_REF] Homsy | Viscous ÿngering in porous media[END_REF][START_REF] Quintard | Stabilità e d'une zone instationnaire de di usion thermique dans un milieu poreux: Analyse linà eaire, modà elisation numà erique[END_REF][START_REF] Quintard | Convection naturelle en milieu poreux: Syst emes non-stationnaires, dà eplacements[END_REF][START_REF] Tan | Stability of miscible displacements in porous media: Rectilinear ow[END_REF], 1992;[START_REF] Yortsos | Dispersion driven instability in miscible displacement in porous media[END_REF] like the frozen time method, the linear amplitude equation method or the energy method that are very brie y presented below.

(i) Frozen time or quasi-static or quasi-steady-state analysis. As mentioned above, the bifurcation theory between steady-state ows is no longer valid when W is a function of space and time. However, assuming that t is a parameter in W (x; t),-this is a strong hypothesis in this theory and is di cult to justify-the base ow concentration ÿeld becomes frozen, and the linear analysis can be formally applied to the disturbances. A lot of sta-bility studies [START_REF] Foster | Stability of homogeneous uid cooled uniformly from above[END_REF]Robinson, 1967;[START_REF] Gresho | The stability of a uid layer subjected to a step change in temperature: Transient vs. frozen time analyses[END_REF][START_REF] Tan | Stability of miscible displacements in porous media: Rectilinear ow[END_REF] have been based on this method since it is easier to implement than the two other methods described below. However, it is of limited usefulness for most transient cases. In particular, the marginal stability, which corresponds to zero time derivative for the disturbance is certainly the worst situation compared to the ÿnite time variations of the base ow. Therefore, it is better to use the conditions for rapid disturbance growth as instability criteria [START_REF] Tan | Stability of miscible displacements in porous media: Rectilinear ow[END_REF]. However, in our case the base state is periodic in time, with complex variations due to the various changes in the boundary conditions. Therefore, it is di cult to estimate the validity conditions for this quasi-steady-state approach.

(ii) Linear amplitude equation method. Linear amplitude equations can be derived from the original governing equations in the case of non-autonomous systems. Solutions of the resulting equations give valuable information about the possible development of perturbations and, in this respect, yields better results than the previous theory. However, it is not suitable for describing well-developed perturbations, which can interact in a non-linear fashion [START_REF] Malher | Stability of a uid layer with time dependent density gradients[END_REF][START_REF] Quintard | Stabilità e d'une zone instationnaire de di usion thermique dans un milieu poreux: Analyse linà eaire, modà elisation numà erique[END_REF].

(iii) Energy method. This method can account for well-developed perturbations [START_REF] Joseph | Stability of uid motions[END_REF][START_REF] Homsy | Global stability of time-dependent ows: Impulsively heated or cooled uid layers[END_REF]. In order to determine the stability conditions, a Liapounov functional is built from the governing equations. This method is interesting since non-linearities present in the system are taken into account. In particular, this method can show how a perturbation can actually vanish when combined with time evolution of the regular solution while it is considered as an unstable one for the two other methods [START_REF] Quintard | Convection naturelle en milieu poreux: Syst emes non-stationnaires, dà eplacements[END_REF]. In other words, this analysis provides critical times at which perturbations may increase or decrease depending on the evolution of the base ow and other parameters [START_REF] Quintard | Instabilità es de zones de di usion thermique instationnaires en milieu poreux[END_REF]. However, it might be di cult to obtain the solution in the case of very complicated problems, i.e., many coupled, non-linear equations. Therefore, the e ect of the non-linearities is often investigated by solving the full set of non-linear equations for various disturbance patterns. This is used to validate the conclusions of the stability analyses. However, it generally requires heavy computations, and cannot be used for a parametric study.

Because of the cyclic nature and the complexity of the system under study in the present work, linear amplitude equations were derived and solved numerically for the stability analysis. To simplify the mathematical treatment of the problem, we assumed that the temperature disturbance has a negligible e ect on the overall stability. This implies that the energy equation will not be considered in the stability analysis, i.e., it will be reduced to Darcy's law combined with the mass conservation. and a can be identiÿed to the wave number as a = 2 : (59)

Initial and boundary conditions

To complete the system of equations, Eqs. ( 46) -( 48), initial and boundary conditions are required and their choice is of capital interest in the stability analysis of non-autonomous systems. In fact, in the case of classical bifurcations between steady solutions, it is su cient to know that perturbations exist. If the system is under unstable conditions, the solution eventually switch from the stable base ow to another solution and the structure of the perturbed state does not often depend on the choice of the initial disturbance. Physically, this means that the source of the instability is irrelevant for the analysis. This statement has to be slightly moderated, however, in the case of sub-critical conditions since the jump from one solution to another may depend in that case upon the disturbance amplitude. This is much more complicated when the system is non-autonomous. The behavior of the system is highly dependent on sources giving rise to the perturbations. For instance, in the case of the stability of transient thermal di usion zones in porous media, it was found that the observed critical times depend on the choice of the initial perturbations [START_REF] Quintard | Convection naturelle en milieu poreux: Syst emes non-stationnaires, dà eplacements[END_REF]. In some cases for which the di usion zone is convected through the porous medium, i.e. a base ow with a non-zero velocity, the perturbation seems to be of limited growth compared to the case of a stagnant zone if we observe the ow over a limited length. This is simply because convection ushes the disturbance away from the domain before it has reached a ÿnite amplitude accessible to observation [START_REF] Quintard | Criteria for the stability of miscible displacements through porous columns[END_REF]. All these aspects are illustrated and discussed below with the presentation of results on the particular problem studied in this paper.

From a mathematical point of view, several di erent choices are possible for the perturbation sources. For instance, the source may be an initial perturbation of ÿnite amplitude. If the system is unstable for the conditions under study, we may observe an increase in the perturbation amplitude. However, and this is especially true for cyclic conditions, the disturbance may be carried away from the domain under consideration by convection, and if there is no source associated with the boundary conditions the disturbances will eventually vanish. Subsequently, there will be virtually no perturbation sources for the following cycles. In a real system, perturbation may originate from sources that do not support the scheme of a ÿeld only perturbed initially. For instance, if the perturbation is initiated by small uctuations of physical properties-like heterogeneities-one can easily imagine that this perturbation is kept all along the cycle and not only at the early stage of the ow. For this reason, a stationary disturbance could be maintained for instance at the entrance, in addition to an initial perturbation. Alternatively, a perturbation could be superimposed to the solution at each time step in a random manner. This last idea was not investigated in the present work and results obtained in the case of an initial perturbation on the one hand and for a stationary disturbance at the boundary on the other hand are presented. This allows the discussion of the stability problem on a general basis keeping in mind that detailed practical applications would require more knowledge of the possible perturbation sources.

Stationary boundary conditions of the perturbed problem were taken similar to the boundary conditions of the regular problem. For instance, Neuman conditions for the regular concentration, Eqs. ( 15), ( 18), ( 21), ( 24), ( 27), ( 30) and ( 36), remain Neuman conditions for the disturbed ÿelds, while Dirichlet conditions for regular concentration and regular pressure, Eqs. ( 11) -( 13), ( 19), ( 22), ( 28), ( 33), ( 37) and ( 39), remain Dirichlet conditions. These Dirichlet conditions were prescribed as follows (B:C:31) Cÿ1 (t; y) = 0 on e (Period I);

where e and s represent the entrance and exit surfaces. The boundary condition in Eq. ( 62) takes into account pressure uctuations occurring at the entrance of the domain. In the case of impervious boundaries, the regular velocity is zero, as expressed by the boundary conditions in Eqs. ( 16), ( 25), ( 31), (34), and this implies that the pressure perturbation is also equal to zero.

In the case of an initial perturbation ÿeld, two different concentration conditions were investigated: an over-oxygenation and an under-oxygenation.

Discretization and numerical procedure

A classical ÿnite volume method was used to solve the regular problem, Eqs. ( 1), ( 3), ( 4), ( 6), ( 7) and ( 9), and the perturbed one, Eqs. ( 46) -( 48), with their corresponding boundary conditions. Integration was performed on the control volume depicted in Fig. 4 where we have also indicated the positions of the unknowns. represented the results for concentration versus the dimensionless position along the porous column axis. These results were obtained during the feeding period after a dimensionless time of 0.033 and 0.133 in the ÿrst cycle. As a comparison, the regular solution was also computed using a 2D numerical model 1 over the same period of time and concentration results obtained at a ÿxed value of y corresponding to the mid-section are also depicted in Fig. 5.

One can observe a small di erence between the two adsorption waves due to the fact that the 2D numerical model takes into account the time to reach adsorption equilibrium, i.e., it uses a local non-equilibrium model, whereas adsorption equilibrium is assumed to be instantaneous in the 1D model. However, the agreement between the two approaches is good, and the regular solution obtained from the 1D model was considered to be accurate enough for the subsequent stability analysis. For the same ow rate, results on oxygen mass fraction and temperature evolution over a complete cycle are represented in Figs. 6 and7 versus dimensionless time and position along the column axis, respectively.

To complete the regular ow analysis, it is also necessary to test the e ect of the cyclic nature of the process on the regular solution. Results on the regular concentration evolutions along the porous column during the early stage of the feeding period over four cycles are depicted in Fig. 8. They indicate that: (i) A cyclic steady state is reached after few cycles (typically 2 or 3) depending on the ow rate. This is in agreement with previous reported results obtained un-1 The 2D numerical model was developed in RAMPANT TM by Fluent Inc.. der isothermal conditions [START_REF] Lu | Simulation of a three-step one-column pressure swing adsorption process[END_REF]. As a consequence, cyclic steady-state conditions are considered as initial conditions in the subsequent stability analysis.

(ii) The porous bed recovery decreases with increasing ow rate. (iii) Steady state is reached more rapidly for low ow rates. In practice, however, gas production rate may require high ow rates despite lower ones are in favor of higher bed recovery and earlier steady state.

Starting from these results, the perturbed equations were solved according to the scheme described above.

Perturbed case

Before presenting our results, it is important to emphasize that the unstable behavior is intrinsic of the system under consideration and is not a pure artifact resulting from the 1D linearization of the system of equations. As an illustration, a 2D numerical solution was obtained on the feeding period using RAMPANT TM with a perturbation for the pressure boundary condition given, as in Eq. ( 62), by Pÿ (t; x = 0; y) = P sin 2 y ;

(66) and a concentration initial perturbation of the form:

Cÿ1 (t = 0; x; y) = Cx(L -x) sin 2 y ; 0 6 x 6 L:

(67)

Pressure results obtained with P =20 Pa, =L and C = 0:05 at t = 0:133 during the feeding period are reported in Fig. 9.

Clearly, the perturbation grows and spreads inside the porous column conÿrming that disturbances may develop inside the medium and justifying a more thorough computation over a complete cycle with the 1D simpliÿed model.

To begin with, perturbed equations were solved with an initial concentration perturbation only, i.e., P = 0, and no disturbances applied on the boundaries. Two di erent types of concentration perturbations were investigated: one corresponding to an over-oxygenation of the column, the other one to an under-oxygenation. Both disturbance ÿelds were of the form given in Eq. ( 67). Computation was performed on the overall cycle and results obtained with C = ±0:05 and = L are reported in Fig. 10 were we have represented the time evolution, over the ÿve periods, of the relative norm of the oxygen concentration, Cÿ1 (t) r; 2 , deÿned by

Cÿ1 (t) r; 2 = Cÿ1 (t) 2 Cÿ1 (t = 0) 2 = L 0 C2 ÿ1 (x; t) dx L 0 C2 ÿ1 (x; t = 0) dx : (68)
It is clear from these results that the perturbation evolution di ers from step to step and depends on the initial perturbation. During purge and delay, the perturbed concentration increases, whatever the initial perturbation. During elution and pressurization steps the perturbed concentration decreases. On the contrary, one can observe a di erent behavior during the feeding period since perturbation increases with an over-oxygenation while it decreases with an under-oxygenation. Although its evolution depends on the velocity inside the column, (see results in Fig. 11 for C = -0:05 and two di erent ow rates), perturbations are always dampened after few steps.

In fact, in this problem, perturbations are carried away from the porous medium by convection while a zero perturbation condition at the inlet tends to be convected into the porous domain since there is no perturbation maintained at the inlet. Therefore, since we are considering a cyclic behavior it is expected that the system is relatively stable. Here the word stable has to be understood in terms of transient analysis. Indeed, several calculations showed that the ow is stable with these boundary conditions and for the physical situation under consideration, whatever the initial concentration perturbation and the wave number. In terms of practical applications, such a situation may seem to be too drastic. In order to release this constraint, it is interesting to investigate the case of a perturbation maintained at the inlet of the column. This kind of condition is physically justiÿed by experimental observations of velocity uctuations at the entrance during feeding, purge and elution steps.

Computations were performed with a perturbed pressure boundary condition of the form given in Eq. ( 62). Results on Cÿ1 (t) r; 2 versus dimensionless time obtained with a ÿxed wavelength, , equal to L, and P ranging from 5 to 30 Pa are represented in Fig. 12. These results clearly highlight the conditional stability of the system. For this wavelength, for instance, the concentration disturbance begins to grow for an amplitude of the pressure perturbation greater than 5 Pa.

However, disturbance amplitude at the entrance is not the only factor in uencing stability since the wavelength is also a key factor conditioning the stability. In Fig. 13, are reported our results for Cÿ1 (t) r; 2 versus dimensionless time, obtained with the inlet pressure disturbance amplitude ÿxed to 20 Pa and wavelengths ranging from 0:13 × L to 10 × L. These results indicate that stability occurs for relatively large wavelengths (roughly larger than 5L in our case), while the time to reach the instability decreases with decreasing wavelengths. For very short wavelengths, the solution diverges very quickly and this directly results from the linear character of the operator.

To synthesize the above results, we have summarized the stability analysis in Fig. 14 providing the stability diagram in terms of amplitude and wavelength of the pressure disturbance at the entrance of the column. This diagram was obtained for a pressure drop of 600 Pa along the column during the feeding period.

As shown in this ÿgure, unstable ow results from large enough disturbance amplitudes at the entrance, i.e., for a su ciently large amount of energy input in the system. The amount of energy required to reach this unstable behavior is an increasing function of the wavelength.

Conclusions

In this work, stability during cyclic gas ow with dispersion and adsorption in a porous column, as encountered during pressure swing adsorption, was investigated. On the basis of a 1D linearized amplitude equation, and from the computed solution of the time-dependent regular ow, a numerical solution was obtained for the concentration perturbation in the case of two di erent kinds of disturbance. For an initial concentration disturbance, it was found that the system is always stable. This result can be explained if one considers the e ect of convection through the porous domain. In fact, since no disturbance is maintained at the entrance, a zero disturbance tends to be convected through the column so that the initial disturbance is ÿnally dampened away. A somewhat di erent behavior was observed if a perturbation is maintained at the inlet for the pressure boundary condition for instance. In this case a conditional stability is obtained on the adsorption front inside the porous medium. It was found that instability occurs for large enough perturbation amplitude at the inlet, and that the amplitude required to destabilize the dispersion=adsorption process increases with increasing wavelengths. Although the in uence of all the physical parameters has not been investigated since it was beyond our scope in the present work, a stability diagram in terms of amplitude and wavelength of the pressure perturbation at the inlet was obtained.

Fig. 1 .

 1 Fig. 1. The three main steps in the adsorption cycle.

Fig. 2 .

 2 Fig. 2. Entrance pressure evolution during a complete cycle.

Fig. 3 .

 3 Fig. 3. Geometry of the macroscopic medium.

(

  C:34) Pÿ (t; y) = 0 on s (Periods I; II; V);

Fig. 6 .

 6 Fig. 6. Regular oxygen mass fraction along the column over a complete cycle.

Fig. 7 .

 7 Fig. 7. Regular temperature along the column over a complete cycle.

Fig. 8 .

 8 Fig. 8. Regular concentration evolution along the porous column at the early stage of feeding for four cycles and two di erent ow rates. (a) V ÿ = 10:854L= ; (b) V ÿ = 2:508L= .

Fig. 9 .

 9 Fig. 9. Simulation with a 2D model: perturbed pressure in the porous medium; t = 0:133 ; V ÿ = 2:508L= .

Fig. 10 .

 10 Fig. 10. Evolution of the relative norm of the perturbed concentration over one cycle resulting from initial concentration perturbations; V ÿ = 2:508 L= .

Fig. 11 .

 11 Fig. 11. Evolution of the relative norm of the perturbed concentration along the column for an initial under-oxygenation and two di erent ow rates.

Fig. 12 .

 12 Fig. 12. Evolution of the relative norm of the perturbed concentration resulting from a perturbed pressure at the entrance. E ect of the perturbation amplitude.

Fig. 13 .

 13 Fig. 13. Evolution of the relative norm of the concentration disturbance resulting from a perturbed pressure at the entrance. E ect of the perturbation wavelength.

Fig. 14 .

 14 Fig. 14. Hydrodynamic stability threshold as a function of the disturbance wavelength and amplitude, for a ÿxed feeding pressure drop.

  in the gas phase, m 2 =s. Here, D * ÿ = D * ÿ I f(C ÿi ) function describing the adsorption isotherm for component i K porous media permeability, m 2 . Here, K = KI K * thermal conductivity of the uid=solid system, W=m K -1 . Here, K * = K * I L column length, m m i mass component i in a mass m of mixture, kg m mass of the gas mixture, kg M ÿi molar mass of component i; kg=mol n i number of molecules of component i contained in a mixture of n m molecules n mnumber of molecules in the gas mixture P C p ) * volumetric heat capacity of the uid=solid system, J=m 3 K -1 ( C p ) volumetric heat capacity of the uid -= ÿ-or solid -= -phase, J=m 3 K -1 complete cycle time, s e ; s entrance and exit surfaces of the bed

Acknowledgements

Financial support from the Air Liquide Company is gratefully acknowledged.

In the following, we use P 0 ÿ ; T; C 0 ÿ1 and V 0 ÿ for the solution to the base ow corresponding to the initial boundary-value problem described by Eqs. ( 1) -( 30). The disturbance ÿelds, Pÿ ; Ṽÿ ; Cÿ , are deÿned as Pÿ = P ÿ -P 0 ÿ ;

(43)

The disturbance linearized equations were obtained by introducing these decompositions into the governing equations, Eqs. ( 1) -( 7). During the linearization process, products of disturbances of the form X ÿ Ỹ ÿ were discarded and the following disturbance equations were obtained:

Darcy's law:

Oxygen mass balance:

Mass balance for the mixture: @ ˜ ÿ @t compression

In the two above equations, a series of approximations were used. In particular, adsorption and compressibility terms were developed on the basis of the following Taylor's expansions:

in which higher second-order terms were neglected. Moreover, in Eq. ( 47), the following inequalities were assumed in order to simplify the dispersive and convective terms:

Physically, this corresponds to the hypothesis of a negligible density disturbance with respect to the concentration and pressure disturbances. This leads to

in Eq. ( 48).

Since the base ow is assumed to be one-dimensional, and due to the structure of the operator, it is possible to decompose the pressure and concentration disturbances under the form Cÿ1 (t; x; y) = C(t; x)g(y);

(55) Pÿ (t; x; y) = P(t; x)g(y):

Once inserted in the linearized equations, this implies that 2 y g + a 2 g = 0;

(57) in which a is a constant. Since the system can be considered as inÿnite in the y-direction, an elementary solution for g is hence, The base ow and disturbed problems are coupled and involve seven unknowns, P 0 ÿ ; Pÿ ; C 0 ÿ1 ; Cÿ1 ; V 0 ÿ ; Ṽÿ and T , which must be computed at the same time step. To do so, convective, dispersive and accumulation terms present in each set of equations were, respectively, discretized according to a speciÿc scheme. Convective (or hyperbolic) terms, like • ( ÿ V ÿ C ÿ1 ) and • (( C p ) ÿ V ÿ T ), were discretized using an explicit upwind ÿrst-order and conservative scheme. The explicit formulation was chosen in order to get a linear form of the problems, and moreover, to provide an e cient way of limiting numerical di usion. Di usive (or elliptic) terms

were discretized with an implicit scheme based on second-order centered ÿnite di erences. Accumulation terms including compressibility and adsorption e ects, @ ÿ =@t and @ C 1 =@t, respectively, were discretized according to the following scheme: @ ÿ @t = @ ÿ @C ÿ1 n @C ÿ1 @t + @ ÿ @P ÿ n @P ÿ @t + @ ÿ @T n @T @t :

(64)

As a result, the discretized form of Eq. ( 1), for instance, reads 55) and ( 56), and these terms were discretized by an explicit scheme.

For each problem, unknowns were computed as the solution of the linear system formed from the set of linearized discrete equations using a classical Gauss elimination method.

Results and discussion

In all our results presented below, length and time are made dimensionless by the column length, L, and complete cycle time, , respectively. Computations were performed using the values of the physical parameters reported in Table 1.

Regular case

To begin with, the 1D base ow was solved using the above-mentioned numerical scheme. In Fig. 5, we have