Mohammed Alloulah
email: alloulah@outlook.com

Airborne Broadband Ultrasonic Tracking: Algorithms and Architectures for Mobile Devices

Keywords:

This paper describes methods for embedded, real-time airborne broadband ultrasonic tracking. The tracker has been built around the assumption of a mobile operation that is deployed ad hoc and upon demand. In order for this to happen, the embedded, real-time operation of sensor nodes has been emphasized. The efficient signalling designs that make way for multiuser, ad hoc tracking deployment are thoroughly characterized, and shown to perform close to their infrastructure-reliant counterparts. System-level parameterization of tracking is also possible subject to application needs. We then demonstrate that real-time Doppler processing in the airborne broadband ultrasonic modality is possible, whereby velocity inference of mobile nodes is facilitated. Building on advancements from underwater acoustics research, a complex Doppler receiver has tough implications on real-time realizations. We study and characterize these requirements utilizing a highlevel synthesis architectural exploration methodology. This has revealed that it is feasible to implement real-time Doppler tracking for airborne broadband ultrasound using modern reconfigurable fabrics.

INTRODUCTION

T HE tracking of mobile users and their devices indoors is of great significance to a number of application domains. Tracking encompasses initial positioning (or localization) followed by continual monitoring of movement thereafter. Numerous examples can be found, from applications such as spatial file sharing in smart environments, computer-assisted living, robotics, all the way to full-blown motion capture. Many technologies have been shown to supply applications with tracking services. These technologies are complementary, and none possesses all desirable tracking attributes for all situations. The preference for the use of a particular tracking technology is often trumped by application-mandated considerations rather than isolated technological merits.

Widely available, pre-existing infrastructures indoors are useful for applications requiring only coarsegrained tracking. The most faithful example of this is perhaps WiFi fingerprinting techniques [START_REF] Bahl | RADAR: an in-building RFbased user location and tracking system[END_REF].

The provision of accurate, fine-grained tracking has been classically achieved by infrastructure-reliant systems. By heavily instrumenting a confined operational space, technologies such as ultrasound, ultrawide band (UWB), optical, magnetic, and hybrid sensing can meet the requirements of all systems [START_REF]InterSense, Inc[END_REF], [START_REF]Ascension Technology Corporation[END_REF], [START_REF]Ubisense Limited[END_REF], [START_REF] Fontana | Commercialization of an ultra wideband precision asset location system[END_REF]; high-end and low-end alike. However, this comes at the expense of high price tag, installation overhead, calibration effort, and end-user expertise.

These factors rule out the cost-effective deployment of these systems for all but the most specialized groups of users e.g. research labs, industrial plants, animation studios, etc.

Ultrasonic trackers have been shown anecdotally to be popular for the provision of tracking services indoors [START_REF] Mccarthy | Accessible Ultrasonic Positioning[END_REF]. The low-rate of ultrasound translates into advantages in terms of accuracy, cost, and ease of installation. However, narrowband ultrasound poses restrictions on the real-world utility of these tracking systems, especially in scenarios pertaining to multiuser operation, high mobility, and rapid deployment in unprepared environments.

The use of broadband ultrasound in-air was first shown to afford a number of system-level enhancements to tracking in [START_REF] Hazas | Broadband ultrasonic location systems for improved indoor positioning[END_REF]. Since, the notion of wideband ultrasonic sensing was further explored in literature [START_REF] Gonzalez | High-Precision Robust Broadband Ultrasonic Location and Orientation Estimation[END_REF], which consolidates the earlier findings in [START_REF] Hazas | A novel broadband ultrasonic location system[END_REF]. However, airborne wideband ultrasonic signals present two primary challenges which have thus far limited their effective use in real-world tracking systems. These challenges fall under two categories: (1) embedded, real-time realization; and (2) proneness to aggressive Doppler distortions as a result of humanscale movements indoors. Consequently, elaborate interactive applications remain difficult to realize using airborne broadband ultrasound (ABU) for location sensing.

This paper seeks to advance the state-of-the-art in ABU sensing on two categoric frontiers of equal importance to real-world trackers: (1) amenability to embedded, real-time realization; and (2) Dopplertolerant tracking for moving tags. In order for this to happen, the paper considers the problem at hand holistically, and jointly addresses design issues across algorithms and architectures, all targeted to mobile devices.

RELATED WORK AND MOTIVATION

Ultrasonic Tracking Systems

Broadband

Since the early implementations of broadband ultrasound for indoor localization [START_REF] Hazas | A novel broadband ultrasonic location system[END_REF], few systems have appeared in the literature describing results in the same band of operation or further exploring signalling aspects. Work by Prieto et al. evaluates the ranging performance of CDMA acoustic signals in both still air and moving air conditions [START_REF] Prieto | Performance evaluation of 3D-LOCUS advanced acoustic LPS[END_REF]. They report subcentimetre error using relatively short (32 bit) Golay codes. However, they utilized more powerful (and bulky) transducers operating in the primarily audible range of 5-25 kHz. This would have allowed for more favourable signal-to-noise ratios than are possible with compact piezofilms operating in the ultrasonic range (above 30 kHz), and it is unclear how much this SNR improvement contributed to their accuracy.

Another system described by Ure ña et al. [START_REF] Ure Ña | Advanced sensorial system for an acoustic LPS[END_REF] replicates previous signalling [START_REF] Hazas | Broadband ultrasonic location systems for improved indoor positioning[END_REF], except that it substitutes Kasami codes for Gold ones. The system relies on the differential time-of-arrival of signals from five transmitters to estimate one receiver's location. Sporadic measurements are shown but no comprehensive characterization of accuracy within the test volume is conducted. The paper also discusses issues of embedded design, considering sequential and parallel implementations of the Kasami code correlator. The authors report the implementation of single-user, parallel correlator with eight-bit resolution which utilized 2,328 slices on a Spartan-IIE architecture. How the receiver embedded processing logic scales as the number of transmitters increases is not discussed.

Álvarez et al. recommend a limit of about 10 ms of Doppler coherence time to govern code length, in order to specifically address the air turbulence effects they observed outdoors [START_REF] Álvarez | High reliability outdoor sonar prototype based on efficient signal coding[END_REF]. For indoor environments turbulence is less aggressive, yet still present due to air conditioning units and building draughts. In the proposed solution for tracking in the presence of Doppler, we advocate that equalization per-chip mitigates against the symbol coherence time problem. The advantage of this method is that there are no hard limits on code length, transmission power, or transducer efficiency. More flexibility in these parameters allows the design of devices with small form factor and portability-crucial for small, battery-powered sensor nodes and wearable tags.

Yet another work by Gonzalez et al. characterizes the ranging and orientation performance of CDMA acoustic signals [START_REF] Gonzalez | High-Precision Robust Broadband Ultrasonic Location and Orientation Estimation[END_REF]. Two spread spectrum (SS) schemes are compared: direct sequence (DSSS) and frequency hopped (FHSS). Both schemes utilized Kasami codes, a 40 kHz carrier frequency, and a bandwidth of 15 kHz, which is considerably less than the bandwidth in our system (40 kHz). The DSSS chip rate was 10 kHz with 200 chips per symbol, while the FHSS used 2.5 kHz and 50 chips. The evaluated ranging performance reported one hundred measurements at four distances. At the maximum of which (1.24 m), the 95 percentile confidence level of the cumulative error was 4.7 mm. For an 8-element uniform circular array, the 3D orientation performance reported was obtained from few discrete locations by averaging 20 estimate at each. This resulted in: 4.5 • yaw, 3 • pitch, and 3.5 • roll accuracies.

Smart Antenna Sensing for ABU Tracking

ABU's ability to natively sense the three physical quantities that constitute the spatio-temporal relationship of a group of co-located users-range, speed, and orientation-presents exciting new opportunities and a set of accompanying challenges. The implications of different ABU sensing topologies on PHY and MAC design are discussed next.

We first remark that in ABU tracking networks, the use of an RF back channel for loose synchronization is assumed. For a group of co-located users, one node will be tasked with regulating sensing through the use of a broadcast-style RF trigger, synonymous to a "tick" in simulation. This seemingly added overhead is a small price to pay for enhancing the longevity of the network and simplifying ABU transceiver requirements. Accomplishing the same task through acoustic messaging will be significantly more expensive both in terms of energy and signalling complexity. It is also assumed that nodes will report their ABU measurements over RF at a much sparser periodicity.

An ABU network consists of a number of nodes, possibly of heterogeneous architectures. A node capable of more functionalities than simple ranging will be designated as basestation (BS). In the following, we further describe the operation of possible topologies, and discuss their differing cross-layer implications.

Topology 1: Uplink

As shown in figure 1a, here ABU tags are transmitters and the BS is a receiver. This topology has an obvious ABU capacity limit, albeit one enhanced by spatial combining. There are three implications topology 1 has on the system-level design.

Range-velocity uncertainty. Though theoretically feasible, multicode acquisition in the dual presence of range and velocity uncertainties is hard and expensive. The relaxed scenario which is multiuser tractable is for ABU messaging to be initiated from a stationary position. One way to reconcile this constraint with usage scenarios is for a tag to request Doppler-tolerant tracking from BS, say using a push-button. Higher-order motion moments. Prevalent as a result of human-scale movements, accurate Dopplertolerant tracking is a challenge in the presence of higher-order motion moments i.e. acceleration and jerk. That is, special signal processing techniques need to be developed in order to realize ABU motion inference resilience irrespective of operational conditions.

Power Control. Power control is important owing to the pronounced distance-dependent channel variations in ABU [START_REF] Alloulah | Airborne Broadband Ultrasound: Wireless Channel and Motion Tracking Algorithms[END_REF]. That is, the physics of propagation across the wide band of acoustic frequencies dictates the need for extensive power control-beyond what various modulation schemes can alone alleviate. For an ad hoc ranging system, power control can only be handled by MAC as to achieve the best possible acoustic equilibrium among transmitters at different distances. Additionally, such loose power control could be implemented by the RF messages that report measurements to the end application. Doppler tracking power control, however, could be achieved to within ±0.2 dB in similar vein to established techniques in commercial systems [START_REF] Holmes | Spread Spectrum Systems for GNSS and Wireless Communications[END_REF].

Topology 2: Downlink

Figure 1b depicts a downlink ABU topology in which tags are receivers and the BS is a transmitter. The tethered BS which has less restrictions on form factor and power affords a drastically different PHY approach. As such, the tethered BS assigns many unique spreading codes, each steered towards a particular spatial sector, sums the resultant beams, and transmits the aggregate beam continuously. The continuous stream resembles that of GPS in that it consists of data frames with acquisition preambles. Consequently, the need for RF synchronization is alleviated altogether. Such topology is infinitely scalable. However, this comes at the expense of equipping receiver tags in topology 2 with binary Doppler-tolerant reception capabilities.

Binary Doppler-Tolerant Reception. This receiver effectively doubles the computational complexity i.e. silicon area. It has been reported in the UWA literature. Dubbed the hypothesis feedback equalizer [START_REF] Stojanovic | Multichannel Detection for Wideband Underwater Acoustic CDMA Communications[END_REF], it entails having two parallel adaptation branches all the time scanning for two [binary] hypotheses ahead of the end of any given symbol in order to mitigate against data transitions.

Topology 3: Hybrid P2P

Another hybrid topology is possible as illustrated in figure 1c. Here the emphasis is on peer-to-peer (P2P) sensing in which a combination of nodes capable of ranging, spacial functions, and Doppler processing is present. This topology relies again on an RF MAC for the successful execution of functionality. Such heterogeneous sensing will arise naturally as a result of the evolution in the network usage over time as nodes opt in and out of the service ad hoc. Thus it is most accommodative of versatile and demanding user experiences. However, this comes at the expense of supporting the reconfigurability of tracking.

Reconfigurability: The ability to adapt tracking to arising needs is a feature that could potentially enhance the user experience. This reconfigurability should at the same time facilitate the full utilization of the ABU channel. That is, the choice of signalling in ABU should lend itself to a multiuser operation that is both computationally tractable and accommodative of different sensing functions; be them ranging, Dopplertolerant tracking, or beamforming for active steering or direction-of-arrival estimation.

Limitations of Ultrasonic Tracking Systems

Broadband

When so many systems in the literature are described by being broadband while having differing bands of operation, performance comparison becomes unreliable. However, three major observations on the stateof-the-art in broadband ultrasonic tracking hold:

1) Unfavourable transducers: Many of the aforementioned systems utilize bulky transmitter and receiver devices which severely limit the usage scenarios of broadband ultrasonic tracking systems. This is especially true for the targeted usercentric applications operating in a mobile fashion, which require light-weight and rugged tags.

That is, the design of signalling techniques for ABU should take into account realistic, imperfect conditions e.g. the low sensitivity of rugged, lightweight, and wideband piezo films. 2) Decreased bandwidth: Due to physical transduction limitations, the above systems utilized only a fraction of the bandwidth originally reported in [START_REF] Hazas | Broadband ultrasonic location systems for improved indoor positioning[END_REF]. Among other things, this influences the capacity of the channel-a property that substantially impacts the perceived quality of service (QoS). This is because both the number of coexisting users (multiple access) and their update rates are increased with the enhancing of capacity as a result of wider bandwidth. 3) Motion-intolerance: All reported systems fail the instance tags begin to move while emitting or receiving the wideband acoustic signals. This inability to perform sensing is due to the pronounced Doppler effect in the ABU band as a result of human-scale movements. This singular limitation severely affects human-centric tracking since stationary sensing can not be guaranteed.

Joint system-level design

Previous approaches do not consider the system-level problem at hand holistically. That is, no attention is paid to how embedded peer-to-peer sensing between mobile devices could be realized. Nor do these works address the problem of extending their choice of signalling to supply say coherent processing in situation when such sensing enriches interactive scenarios. For example extending FHSS to coherent processing for velocity estimation is, while being theoretically feasible, difficult [14, p. 57]. Hence approaches that touch on one aspect of a tracking system in isolation of equally important system-level issues are bound not to take off. Another example on this point is the use of super-resolution array processing methods that give degree-level accuracy. This level of accuracy is irrelevant when a realistic system performing peer-topeer sensing ought to support multiuser beamforming with potentially a modest but sizeable number of users vying for the service. With a cross-layer approach, the new notion we advocate for is involved processing for the novel ABU modality. The aim is to support a new class of emerging, infrastructureless systems that place particular emphasis on uninstrumented and unprepared environments through self-contained P2P sensing. However, the area-power-cost requirements of high-complexity wireless signal processing are best addressed through joint optimization, wherein algorithms and architectures are cross examined in order to arrive at a tradeoff particular to a set of design constraints [START_REF] Austin | Mobile supercomputers[END_REF], [START_REF] Arvind | Highlevel synthesis: an essential ingredient for designing complex ASICs[END_REF]. For high-complexity algorithms, joint optimization also calls for an agile electronic system-level (ESL) design methodology. Under such a methodology, algorithmic-architectural design space co-exploration is steered judiciously by the informed user while being performed automatically with guaranteed quality of results [START_REF] Cong | Customizable domain-specific computing[END_REF], [START_REF] Guo | Rapid industrial prototyping and SoC design of 3G/4G wireless systems using an HLS methodology[END_REF].

This research reports on advances for ABU PHY designs targeted at multiuser ranging and Dopplertolerant tracking. These designs are proposed while keeping a system-level view at gaze for a rapidly deployable ABU node topology in unprepared environments. In order to achieve the objective of rapid deployment, we outline our signalling approaches alongside their requisite digital architectures, demonstrating how to realize an end-to-end system.

ABU TRACKING

The proposed ABU tracker is comprised of a transmitter with configurable parameters and two modes of receiver detection; static and Doppler-tolerant.

Transmitter

The transmitter is capable of changing its spreading rate as to increase the ABU channel spectral occupancy. Also depending on the receiver mode, it either sends a burst of one code for time-of-arrival (TOA) estimation or a continuous packet transmission for velocity estimation through the Doppler effect. A block diagram of the transmitter is shown in figure 3. The digital transmitter consists of a Gold code generation unit, a numerically-controlled oscillator (NCO), and a DSSS modulator. A mixed-signal DAC then follows whose output is in turn smoothed with an analogue reconstruction filter. Various signal scaling, conditioning, and amplification are applied before finally driving the piezo film.

In order to highlight the signalling adaptability of the transmitter, two chip rate transmission schemes are contrasted next. The nominal 20 kHz (or kChips/s) chip rate was doubled and a raised-cosine filter with a roll-off factor of 0.15 was used to limit the bandpass bandwidth so as to operate within the characteristics of the combined transmitter-receiver ABU wireless channel. The remaining system parameters were: 200 kHz sampling frequency for both the transmitter and the receiver, and 10x & 5x oversampling for 20 kHz & 40 kHz chip rates, respectively. Under this setup, the PSDs of the transmitter's excitation signals are illustrated in figures 2a and 2d. It is clear that the doubling of the chip rate spreads the energy of the DSSS ABU signal more around the carrier (50 kHz).

In relation to the ranging performance, a transmitter-receiver pair is positioned on-axis at a distance of 2.741 m. The time-of-flight of a codelength burst is measured under the two chip rates. A one-second continuous transmission is sampled at the receiver for the evaluation of the PSD for both chip rates as well. Figure 2 depicts the correlation and PSD obtained from the two configurations. The calculated distances from the correlation peaks of the 729 m, respectively. The ranging errors for both cases are on the order of 1 cm. Inspecting the corresponding correlations of figures 2b & 2e reveals that the double rate correlation gain is reduced with respect to the single rate. This reduction is around 6.5 dB. However, so does the noise of the correlation-the faster rate seems to have "dithered" the noise better.

As expected from Parseval's duality principle, the 40 kHz PSD in figure 2f is also reduced similarly with respect to the 20 kHz PSD of figure 2c. It also can be seen that the double rate PSD is more spread around the 50 kHz carrier but remains buried in the noise floor of the receiver outside the intended range due to stricter roll-off factor control.

Receiver

There are two detection modes that an ABU receiver needs to support:

• Static: In this case, nodes are stationary most of the time and rely on TOA inference after a broadcast-style RF trigger has been sent to initiate the pseudoranging1 interval for a group of colocated users. of the continuous wave will inevitably interfere with the direct path and degrade the measurement setup. Both issues are overcome in the ABU method proposed in this thesis. First, by making a "standing" spreading code equate to the maximum expected distance, code timing readily achieves distance disambiguation, of course provided that transmission and acquisition commence at the same time-The RF broadcast-style trigger message will guarantee this. Second, in direct sequence spread spectrum (DSSS) signalling, multipath arrivals at the receiver are resolvable, alleviating the problem in narrowband waves.

Static

The ranging performance of a static, embedded ABU receiver is illustrated in figure 4. The cumulative distribution of ranging errors was generated using a dataset gathered from real ABU deployment through cosimulation. As such, results are representative of real-world, resource-constrained performance.

The experimental conditions of the measurement setup are summarized here. Gold codes of length 511 were used to BPSK modulate a 50 kHz carrier, at a chip rate of 20 kHz. The signal from each receiver's analogue stage was sampled at just over 200 kHz, with twelve bits of resolution per sample.

A primary transmitter node was used to emit one hundred ranging messages, across sixteen horizontal positions (shown as diamonds in Figure 4a), at each of three different heights (10, 50 and 100 cm above the floor). To emulate the presence of co-located users, four other transmitter nodes (triangles in Figure 4a), emitted their ranging messages simultaneous to the node placed at the test positions. All transmitters emit ranging signals at the same acoustic power. Four receiver nodes (stars in Figure 4a) were mounted on the ceiling, 220 cm above the floor.

Figure 4b illustrates the multiuser performance of the co-simulated dataset. For each of the four ceiling receivers, the one hundred ranging events were processed at the forty-eight transmitter locations-a total of 19,200 measurements. The core's overall detection rate for ranging signals from the primary transmitter node is 65.4%. 2 The undetected ranging signal occurrences are commensurate with two effects previously observed with this deployment: (1) the directivity of the transducers [7, figs. 7 & 11]; and (2) transmitter near-far effects [7, sec. 4.3]. When the despreader core successfully detects events, its ranging error is better than three centimetres (2.92 cm) at the ninety-fifth percentile confidence level.

Doppler-tolerant

Tracker configurations. There are three configurations that demonstrated ability to track the Doppler effect in ABU: the phase-locked loop (PLL), decision-feedback equalizer (DFE) with embedded PLL, and DFE with linear interpolator (LI). The following will quickly review these configurations.

(1) PLL: In and by itself, the PLL does not constitute a reliable tracker configuration. However, it is instructive to study the architectural implications of the PLLalone tracker since the DFE-PLL configuration does result in motion-correlated measurements. As such, cross examining the PLL and DFE-PLL tracker configurations will build an understanding on how the architectural requirements change with the inclusion of the feedforward filter (FF) which completes the implicit DFE-PLL tracker.

(2) DFE-PLL: The DFE-PLL tracker configuration provides a means for directly estimating the instantaneous frequency (IF) resulting from motion stresses in ABU. It could be of particular merit in certain situations wherein the general pattern of motion is 2. Ranging success rates higher than the reported 65% would be easily achievable using dynamic thresholding and enhanced peak detection methods [ sought after, rather than accurately measuring velocity. That is, this configuration trades accuracy for robustness. Interestingly though, the DFE-PLL tracker is only valid for a chip oversampling rate of N s = 4. This is a purely algorithmic finding which was arrived at after thorough empirical evaluation [START_REF] Alloulah | Real-Time Tracking for Airborne Broadband Ultrasound[END_REF].

(3) DFE-LI: The DFE-LI configuration is the most accurate velocity tracker. This, however, comes at the expense of sensitivity to high order motion moments (e.g. acceleration and jerk), which are prevalent in ABU as a result of human-scale movements indoors.

Signal Design and Experimental Setup. A Gold code sequence of length 511 was BPSK-modulated onto a 50 kHz carrier at a 20 kHz chip rate. Chips were pulse-shaped using a raised-cosine filter with roll-off factor of 1. The sampling rate utilized was 160 kHz.

One dataset will be used for the characterization of the two algorithmic configurations. This dataset has typical maximum speed (about 0.8 m/sec) and contains high-order motion moments. In order to facilitate characterization, we utilize a vision-based tracking system for the generation of real-time position and speed groundtruth. A transmitter is mounted on a Lego Mindstorms robot moving forward and backward on a track while facing a stationary ABU receiver, at a distance of 2.9 m. The robot is tagged with a fiducial marker [23]. A PC-based application monitors the movement of the robot using the computer vision-based tracking. Upon motion, the application triggers an FPGA-based ABU transmitter and simultaneously commences data acquisition. The alldigital transmitter has an independent oscillator to the PC clock. The experimental setup used for data collection is shown in figure 6.

Performance. The acquired signal was down-mixed to baseband where all processing takes place. Figure ?? illustrates the performance of the DFE-LI tracking configuration. As evident from the first plot, the channel is stable and code timing is maintained throughout tracking. The second plot depicts of the evolution of the linear interpolation factor in time and the chip estimate scatter. The sensitivity of the DFE-LI tracker to acceleration and jerk was avoided initially by a priori knowledge of the Doppler stimulus. As such, the instance at which LI is activated was delayed. As shown, DFE-LI tracks Doppler very accurately during the first and most of the second motion legs. The velocity estimate for the first motion segment was measured to be 0.8047 m/s, which is virtually identical to groundtruth. However, following the decelerating motion moments in the second leg, LI diverges away from rest condition owing to high second-order motion stresses. The third plot examines this stress resulting from second-order motion moments. The IA estimator shows distinct sustained acceleration spikes in three places: at the beginning of the Doppler stimulus around 0.5 sec; around the momentary pause of the robot at roughly 3.5 sec between leg 2 & leg 3; and at the end of the test case slightly before 7 sec. These instances are in line with the acceleration groundtruth of the Doppler stimulus.

Similarly, the performance of the DFE-PLL tracking configuration is illustrated in figure ??. Justifiably, the crude DFE-PLL tracker is noisier even though code timing stays in-lock. In the second plot the PLL phase evolution is shown alongside the chip estimate scatter, which is considerably more spread in comparison to that of the DFE-LI configuration. The third plot shows the scaled unwrapped phase of the chip estimate overlaid on the vision tracking groundtruth. The scaled instantaneous phase is evidently correlated with motion in this competitive mode, which provides a crude means for estimating velocity. This method for inferring the motion pattern could be useful in certain applications not requiring precise velocity measurements.

ARCHITECTURES

Transmitter

The transmitter architecture is modelled in both floating-point and fixed-point arithmetic. The objective is to demonstrate that the fidelity of the proposed embedded realization is commensurate with what the airborne broadband ultrasonic channel has to offer.

The stimulus is baseband digital Gold code signal with infinite quantization precision. Using the fixedpoint model of the transmitter, table 1 shows the obtained SNR at two chip rates. Varying the chip rate has a mild effect on the SNR of the transmitter's excitation signal. Roughly a drop of 4.2 dB accompanies the doubling of the chip rate. However, the two SNRs are significantly higher that that of the aggregate ABU channel which stands at approximately 24 dB, at one meter range [21, p. 76]. The signal at baseband is binary, which results in efficient digital structures and similarly simplifies arithmetic operations. The silicon area of the transmitter prototype for both parameter sets analyzed earlier is given table 2. Evidently, the transmitter is compact and lends itself well to low cost, low power realizations for handheld devices.

Static receiver

An efficient core that makes acoustic simultaneous multiple access possible is shown in figure 7. The proposed core consists of the following units: a complex circular buffer, user code generation unit, a matchedfilter bank, and a control state machine. The matchedfilter bank has a number of complex fingers. The complex fingers are aimed at parallel multiuser reception as opposed to multipath. Each finger is assigned two codes in an interleaved fashion to maximize resource sharing. The Gold code is generated on-the-fly and stepped through the operation, rather than storing the long code sequences in memory. To demonstrate the design, an eight-user ultrasonic despreader was implemented on a Xilinx Virtex-4 XC4FX12 device using RTL description. This was carried out using Xilinx's 10.3 ISE tool suite. The core operates on a 100 MHz system clock with data quantised to 16.15 signed fractional representation, and 8 guard bits. Signalling parameters utilised are a 20 kHz chipping rate, 511 bit Gold code (which leads to about 8.7 m maximum node separation in Dopplertolerant mode), and four-times chip oversampling rate. Under these parameters, four fingers are required, since a maximum of two interleaved users can be accommodated for in a single complex finger. Table 3 shows the final place-and-route metrics of the fourfinger, eight-user core.

ULTRASONIC DESPEADER

Doppler-tolerant receiver-Co-exploration

This section describes the experience of utilizing an HLS tool for C++ to RTL generation. The following presents hardware design issues that are specific to the algorithmic workload.

Chip Oversampling

In summing up the spirit of the proposed chip-oriented architecture, spreading code vectors' worth of oversampled entities are fed to the custom computational datapath. These entities can be an oversampled chip, a reconstructed oversampled chip, or a fractional channel coefficient. The algorithm consists largely of vector arithmetics at the DSSS length (L × N s ) and convolution loops at the code length (L) applied to fractional entities. With 20 kHz chip rate and a nominal digital system clock of 100 MHz, the maximum clock budget per adaptation step is 5000 cycles. While this is quite long in digital terms, the number of loops that have to be performed on the DSSS-long vectors can easily exceed a chip duration. This is because a large proportion of adaptation is purely sequential in nature entailing vector data dependencies. Nonetheless, the one commonality among all computations is that operations are replicated on oversampled entities. Thus a form of data-level parallelism underlying processing is readily achievable by the concurrent manipulation of an N s -wide stream. To this end, as will be further explained in subsequent sections, all loops are partially unrolled by N s while memories are partitioned cyclically by the same factor as to be able to service the N s -wide computational datapathhence the chip-oriented architecture.

Parallelism

The algorithm exhibits a strong form of data-level parallelism determined by the chip oversampling rate. However, in order to tap into this underlying algorithmic uniformity, a few other considerations mandated by the processing workload have to be met, such as expression-balancing some timing-critical operations. This rules out the possibility of using a generic SIMD DSP for realizing a computationally-intensive algorithm of such specific nature, while having equally stringent application requirements on form factor and power consumption. Task (functional) parallelism remains a challenge for HLS, mostly because of the inherent sequential semantics of high-level languages [START_REF] Papakonstantinou | Multilevel granularity parallelism synthesis on FPGAs[END_REF]. In the context of the co-exploration reported in this paper, latency is a critical requirement for the overall adaptation, as is often the case for real-time transceiver designs. Thus all functions are hardware-inlined in order to remove function hierarchy and eliminate calling overhead. Consequently, the functional parallelism problem for HLS is reduced to scheduling, which can be inferred by data dependencies analysis.

Instruction parallelism is automatically handled by HLS. Since the early introduction of HLS as a viable design methodology, more credible automation can now extract instruction parallelisms seamlessly, and can be further aided by the judicious intervention of an expert user through code transformations [START_REF] Buyukkurt | Impact of high-level transformations within the ROCCC framework[END_REF]. This allows for added parallelism subject to broader design space investigation for various trade-offs (such as area, speed, etc) by means of which the final solution is steered.

For the real-time adaptive application at hand, boosted levels of instruction parallelism have been further extracted by pipelining the numerous loops that constitute the algorithm. In addition, the original software code underwent a major refinement with respect to instruction parallelism. This involved visualizing a classic pipeline in mind and rewriting the untimed, sequential C++ code as to reflect a pipelined operation with prefetch, execute, writeback if applicable, and update. This is important so that the pipeline initiation interval (II) remains at one even though the pipeline depth itself could be as high as 30 in certain loops, hence resulting in no pipeline stalls. In such a case, after a small initial ramp-up latency (e.g. 30 clock cycles), the throughput of a loop would be maximized at one operation per clock cycle.

Another remark pertaining to parallelism is the following. In the object-oriented design, global parameters were loaded in variables local to loops or functions prior to proceeding with computations. This is necessary for maximum throughput in these rather long loops. Failure to do so results in long pipeline stalls since the tool effectively translates every occurrence of a global parameter into a function call, which greatly degrades operation scheduling.

Memory

Partitioning. Exploiting the underlying data-level parallelism necessitates a multiport memory architecture as to increase data throughput. In order to read out a complete oversampled chip from the signal vectors in one clock cycle, memories were cyclically partitioned by N s . This is done without introducing any modifications to the code by applying high-level directives. Cyclic partitioning of a signal vector results in a number of equivalent smaller vectors whose entries are the parent's entries interleaved across the now increased memory ports by the required factor. This uniform partitioning is applied to all DSSS vectors; circular and linear. Furthermore, both the Gold code sequence and its estimate are similarly partitioned. This is done for a different reason as will be discussed in the following section. Finally, since the feedforward filter is targeted towards chips, it is much shorter and as such does not require partitioning.

Consistent addressing. It is required that the memory access pattern for various vectors be consistent throughout the algorithm. This entailed refining the original code such that addressing these memories becomes unified as either high-to-low or low-to-high. This aids streaming of data from one function/loop to another and is of particular importance when vector data dependencies occur. Abiding by this consistent access pattern ensures that a datum entry can be immediately streamed to where it is next required as opposed to having to wait on the whole vector operation.

Contention.

(1) access stall: Various functions combine to make a heavy use of the two circular signal vectors; namely, the signal and average signal vectors. Aggravated by vector data dependencies at times, this gives rise to access bottlenecks due to multiplexing and the limited memory ports available coupled with operation rescheduling and/or retiming, which the tool also automatically performs part of a global optimization problem. These access bottlenecks resulted in pipeline initiation interval (II) of value 2 for two loops in the adaptive algorithm; effectively, doubling their latencies. Initially, this was mitigated against by fully partitioning the two circular signal vectors into registers. In the fully partitioned architecture, the circular vectors have as many ports as entries. This is obviously of tremendous generic fabric cost not to mention the very long decode logic delays on these rather deep DSSS-length vectors. Synthesis had proceeded successfully initially. It was not until placing-and-routing the design that failure to meet timing objectives under the fully partitioning highlevel directive surfaced. Therefore, memory access bottlenecks and subsequently pipeline stalls on the two circular signal vectors had to be tolerated for these two loops and were deemed unavoidable as opposed to going fully partitioned. It is interesting to note that this is only applicable to the N s = 4 casethe N s = 2 case did not result in access stalls on these circular vectors.

In order to compensate for these bottlenecks, parallel gain has to be extracted from somewhere else. Bearing in mind the extremely low hardware cost of the last loop despreading_loop-i.e. binary coefficients to put it in signal processing terms-it was decided to partially unroll the loop by the same factor that occurs throughout all loop unrolling directives; N s .

(2) writeback dependencies: The sequence of computations involving the channel estimate vector amounts to four loops; namely, updating, finding maximum coefficient, truncating if applicable, and calculating power. In order to achieve maximum throughput on this subset of algorithmic operations, the following techniques were utilized resulting in modifications to the original code. Firstly, a temporary channel estimate vector was allocated and passed to a first of two channel-related functions which manually merges updating and searching for the maximum coefficient. The updated channel estimate was written to the temporary vector as to mitigate against writeback dependencies that would have otherwise arisen in the pipelined loop. Expression balancing is also applied on the update loop. Simultaneously on the fly, a fast complex magnitude estimator [START_REF] Allie | High-Speed Square Root Algorithms[END_REF] is used for the max search. Listing 1 demonstrates the expressionbalanced, pipelined loop. The second loop then takes on the tasks of truncating the temporary channel vector into the original channel vector while calculating the channel power at the same time.

Listing 1 illustrates the various concepts discussed thus far. First, the channel and interference-free signal are fetched from memory into pipleline registers. Second, multiplications are performed and stored in intermediate registers. Third, the channel is updated and the result is written back to the temporary channel vector in order to avoid writeback dependencies. Fourth, at the same time, the magnitude of the current channel coefficient is computed and the result is written back to the channel magnitude vector. Finally, the current magnitude is tested against the maximum coefficient, which is updated correspondingly if applicable.

(3) feedforward filter adaptation: The complex feedforward filter is implemented as an FIR filter with a RAM-based circular buffer. This is more economical than a register file realization especially given the real and imaginary parts of the signal and coefficient buffers; which means quadrupling the generic fabric utilization for such a realization. As a result of the relatively short delayline of the N s -order filter, adaptation was implemented efficiently incurring only twice the filter's length in clock cycles latency. The two delaylines are loaded into local arrays reversing the signal vector. Indexing the delaylines has to be done in a consistent manner with the actual filtering subroutine as elaborated on above. The local arrays are then used with the adaptation factor to update the original coefficient vector while always pipelining and expression balancing for maximum throughput.

Loops

Various loop-related arrangements have already been alluded to in prior discussion. Nonetheless, for the sake of a cohesive treatment of hardware concepts, the transformations pertaining to loops will be enumerated again under this heading. First, channelrelated loops are manually merged due to the reasons discussed above. Apart from this manual step, subsequent operations are performed automatically through the use of high-level directives. Second, unrolling is applied on virtually all loops occurring in the algorithm. Nested loops are fully unrolled. The remaining loops in which intensive processing takes place are partially unrolled by N s . Third, pipelining directives are utilized on the partially unrolled loops for maximal parallelism.

Co-exploration results

Having described the concepts surrounding the algorithm implementation, the results will be supplied next. The results are based on HLS synthesis by Xilinx's AutoESL 3 version 2010.a.2. The device targeted is the Xilinx Virtex5 XC5VSX94T, package FF1136, speed grade -1. The basic synthesis commands are clock period of 10 ns with 1.25 ns uncertainty. The back-end tools are the Xilinx ISE 12.2 suite.

The IO used in the RTL implementation settings are an input complex chip of type RAM and an output flag of type WIRE connected in a slave bus setup. The Figure 8 shows the characterization of area and throughput as reported from synthesis results. Placeand-route (PAR) was used to sanity-check these results. However, PAR results are not be reported. Besides the long completion overhead, this is also because the downstream tools will prune and optimize away all unnecessary logic that is not tied to the final output assignment. The output assignment is currently a single boolean flag indicating [continual] acquisition. A production-level implementation has to communicate out more algorithmic primitives such as chip estimates for estimating velocity, timing for range 4 , channel status for fidelity, etc. Note that when requiring more IO assignments, no additional latency will be incurred since external communication can be overlapped with computations [START_REF] Cong | FPGA-Based Hardware Acceleration of Lithographic Aerial Image Simulation[END_REF] using nonblocking reads/writes [START_REF] Fingeroff | High-Level Synthesis Blue Book[END_REF]. It is also worth nothing that in the context of automated architectural exploration, prior work in literature oftentimes reports on execution cycles and area estimates from HLS only in avoiding the lengthy overhead of synthesis and PAR [START_REF] Ziegler | Evaluating heuristics in automatically mapping multi-loop applications to FPGAs[END_REF], [START_REF] So | Using estimates from behavioral synthesis tools in compiler-directed design space exploration[END_REF].

The rapid co-exploration allowed us to gain indepth insights into the effect the essential algorithmic parameters have on the chip-oriented architecture. Perhaps most interesting of all was the almost noneffect halving the spreading code (L) had on BRAM usage. This is due to the sequential nature of most of the operations in the algorithm. In the proposed architecture, data is streamed from the global data plane to where it is needed. This streaming is achieved by high throughput FIFOs which are implemented using BRAMs. The number of required FIFOs is algorithmdependent as opposed to code length-dependent. As expected, the major effect of halving the code length Fig. 9: CDFG of tracking configurations on the chip-oriented architecture is seen in the doubling of the computational throughput. On the other hand, halving the chip oversampling rate results in mild throughput enhancement. In order to reap area savings as we decrease code length, modifications to the high-level directives of the chip-oriented architecture have to be made: cyclic memory partitioning and loop unrolling can be relaxed by a factor of two, in order to balance attainable throughput with the requisite area.

Figure 9 shows the control data flow graph (CDFG) for all tracker configurations. The constraining effect of the chip-oriented directives is most evident in the PLL and DFE-PLL configurations at N s = 4-a horizontal increase in the scope of concurrency can be immediately spotted. The DFE-LI tracker, notably, has very favourable data dependencies with two distinct execution branches of minimal coupling towards the latter part of adaptation.

MODIFICATIONS NEEDED FOR A PRO-DUCTION SYSTEM

The set of designs reported in this paper should be viewed as a comprehensive feasibility study on the practical deployability of peer-to-peer, ad hoc ABU sensing. In order to realize a production-ready system, a number of further enhancements and capabilities would be necessary to engineer.

The co-exploration of Doppler tracker configurations can be further optimized in terms of fixed precision datatypes. The Doppler effect is non-stationary and a manual fixed-point precision optimization has to be pursued iteratively. This is formally referred to as simulation driven analysis in [START_REF] Bouganis | Synthesis of DSP Algorithms from Infinite Precision Specifications[END_REF]. Thus, the numerical optimization is laborious and should be deferred until a genuine commitment for building a real-time ABU motion tracker is reached. Note that the numerical design optimization in NASA's Apollo mission consumed around 30% of total development resources [START_REF] Kreide | Computation: Aerospace Computers in Aircraft, Missiles and Spacecraft[END_REF].

Analog circuitries that facilitate power control would be needed for both static and Doppler-tolerant tracking as discussed earlier. In this paper, no attempt was made at realizing any form of power control. The keen reader is referred to [START_REF] Alloulah | Airborne Broadband Ultrasound: Wireless Channel and Motion Tracking Algorithms[END_REF] for a thorough analysis on the distance-dependent variations in the ABU band and a rationalization for ABU power control.

The reconfiguration of tracking is to be supported by a mixture of algorithmic parametrization and hardware reconfigurability. Further significant developments are needed to realize a configurable ABU tracking system that can dynamically cater for changing end-user needs. Customizable internet devices are part of an increasing trend to leverage finergrained hardware reconfiguration in support of network adaptability to dynamic conditions.

CONCLUSION

This paper has developed the novel modality of airborne broadband ultrasound (ABU) in support of ad hoc, real-time, and mobile tracking applications.

Fig. 1 :

 1 Fig. 1: Smart antenna ABU networks

Fig. 2 :Fig. 3 :

 23 Fig. 2: Parseval view of the ABU DSSS system at two chip rates: 20 kHz & 40 kHz.

Fig. 4 :

 4 Fig.4: Ranging performance for a static, embedded ABU receiver.

Fig. 5 :Fig. 6 :

 56 Fig. 5: Tracking performances: upper row DFE-LI and lower row DFE-PLL.

Fig. 7 :

 7 Fig. 7: Multiuser despreader block diagram

3 .Fig. 8 :

 38 Fig. 8: HLS metrics comparison between PLL, DFE-PLL, and DFE-LI

4 .

 4 at least once initially (a) PLL| Ns =2 (b) PLL| Ns=4 (c) DFE-PLL| Ns =4 (d) DFE-LI| Ns=2

TABLE 1 :

 1 Transmitter SNR

	Chip Rate	Roll-off (β)	Input SNR	Output SNR
	(kChips/s)			(dB)
	20	1.0	∞	67.64
	40	0.15	∞	63.47

TABLE 2 :

 2 Tx FPGA area requirements a

	Resource	Used	Utilization %
	DSP48s	1	3
	BRAMs	5	13
	Slices	308	5

a Device: Xilinx Virtex-4 XC4FX12

TABLE 3 :

 3 FPGA resource requirements a

	Resource	Used	Utilization %
	Slices	632	11
	LUTs	699	6
	BRAMs	4	11
	Max frequency		103.67 MHz

a Device: Xilinx Virtex-4 XC4FX12

This term is borrowed from the GPS world and is intended to reflect the uncertainty of the measurement process.


Eduardo Casas
email: eduardo.casas@unican.es

Second Order Conditions for L 2 Local Optimality in PDE Control ⋆

Keywords: optimal control of partial differential equations, semilinear partial differential equations, second order optimality conditions, bang-bang controls

In the second order analysis of infinite dimension optimization problems, we have to deal with the so-called two-norm discrepancy. As a consequence of this fact, the second order optimality conditions usually imply local optimality in the L ∞ sense. However, we have observed that the L 2 local optimality can be proved for many control problems of partial differential equations. This can be deduced from the standard second order conditions. To this end, we make some quite realistic assumptions on the second derivative of the cost functional. These assumptions do not hold if the control does not appear explicitly in the cost functional. In this case, the optimal control is usually of bang-bang type. For this type of problems we also formulate some new second order optimality conditions that lead to the strict L 2 local optimality of the bang-bang controls.

Introduction

This paper is split into three parts. In the first part, we consider the following infinite dimensional abstract optimization problem. Let U ∞ and U 2 be Banach and Hilbert spaces, respectively, endowed with the norms • ∞ and • 2 . We assume that U ∞ ⊂ U 2 with continuous embedding; in particular, the choice U ∞ = U 2 is possible. A nonempty convex subset K ⊂ U ∞ is given, and A ⊂ U ∞ is an open set covering K. Moreover, an objective function J : A -→ R is given. We consider the abstract optimization problem (P) min u∈K J(u),

where we assume that J is of class C 2 with respect to the norm • ∞ . In the next section, we will impose some other assumptions on J so that the first order optimality conditions and the inequality J ′′ (ū)v 2 > 0 for every v ∈ C ū \ {0} imply that ū is a strict local minimum of (P) in the U 2 sense. Here C ū denotes the usual cone of critical directions that we will define later. This result is new in the sense that the classical theory claims the local optimality only in the U ∞ sense due to the non-differentiability of J in with respect to • 2 . Moreover, a stronger inequality J ′′ (ū)v 2 ≥ δ v 2 2 is usually required. In the second part of the paper, contained in §3, we prove that the abstract assumptions are fulfilled by a typical Neumann control problem. The method used for this control problem can be extended in an easy form to many other control problems associated with elliptic or parabolic equations; see [START_REF]Second order analysis for optimal control problems: Improving results expected from abstract theory[END_REF]. Finally, the last part of the paper is considered in §4. There, we analyze the case of bang-bang control problems, which do not satisfy the assumptions of §2. For these problems we also give some second order conditions leading to the strict L 2 local optimality of the controls.

An abstract optimization problem in Banach spaces

The results presented in this section were obtained in collaboration with Fredi Tröltzsch. The reader is referred to [START_REF]Second order analysis for optimal control problems: Improving results expected from abstract theory[END_REF] for the proofs and details.

In this section, we study the abstract optimization problem (P) formulated in the introduction. Besides the hypotheses established in §1 on U 2 and U ∞ , we require the following assumptions on (P). (A1) The functional J : A -→ R is of class C 2 . Furthermore, for every u ∈ K there exist continuous extensions

J ′ (u) ∈ L(U 2 , R) and J ′′ (u) ∈ B(U 2 , R), (2.1)
where L(U 2 , R) and B(U 2 , R) denote the spaces of continuous linear and bilinear forms on U 2 , respectively. (A2) For any sequence

{(u k , v k )} ∞ k=1 ⊂ K × U 2 with u k -ū 2 → 0 and v k ⇀ v weakly in U 2 , the conditions J ′ (ū)v = lim k→∞ J ′ (u k )v k , (2.2)
J ′′ (ū)v 2 ≤ lim inf k→∞ J ′′ (u k )v 2 k , (2.3
)

if v = 0, then Λ lim inf k→∞ v k 2 2 ≤ lim inf k→∞ J ′′ (u k )v 2 k , (2.4)
hold for some Λ > 0.

The reader might have the impression that Assumptions (A1) and (A2), mainly (A2), are too strong. However, we will see in the next sections that they are fulfilled by many optimal control problems.

Associated with ū, we define the sets

S ū = {v ∈ U ∞ : v = λ(u -ū) for some λ > 0 and u ∈ K} , C ū = cl 2 (S ū) ∩ {v ∈ U 2 : J ′ (ū)v = 0} D ū = {v ∈ S ū : J ′ (ū)v = 0}, (2.5)
where cl 2 (S ū) denotes the closure of S ū in U 2 . The set S ū is called the cone of feasible directions and C ū is said to be the critical cone. It is obvious that cl 2 (D ū) ⊂ C ū. However, the equality can fail. In fact, this equality is a regularity condition equivalent to the notion of polyhedricity of K; see [START_REF] Bonnans | Second-order analysis for control constrained optimal control problems of semilinear elliptic systems[END_REF] or [1, §3.2]. This property is enjoyed by control problems with pointwise control constraints. Now, we formulate the necessary first and second order optimality conditions. The second order conditions hold under the mentioned regularity assumption; we refer to [1, §3.2] or [START_REF] Casas | A general theorem on error estimates with application to elliptic optimal control problems[END_REF] for the proof.

Theorem 2.1. Assume that (A1) holds and let ū be a local solution of (P) in U ∞ , then J ′ (ū)(u -ū) ≥ 0 ∀u ∈ K. Moreover, if the regularity condition

C ū = cl 2 (D ū) is satisfied, then J ′′ (ū)v 2 ≥ 0 holds for all v ∈ C ū.
Now, we state our result about sufficient sufficient second order optimality conditions. As the reader may check, the gap between the necessary and sufficient second order conditions is minimal, the same as in finite dimension.

Theorem 2.2. Suppose that assumptions (A1) and (A2) hold. Let ū ∈ K satisfy the first order optimality condition as formulated in Theorem 2.1, and

J ′′ (ū)v 2 > 0 ∀v ∈ C ū \ {0}.
(2.6)

Then, there exist ε > 0 and δ > 0 such that

J(ū) + δ 2 u -ū 2 2 ≤ J(u) ∀u ∈ K ∩ B 2 (ū; ε). (2.7)
Above B 2 (ū; ε) denotes the ball of U 2 with center at ū and radius ε. This theorem can be proved arguing by contradiction. To this end, we assume that for any positive integer k there exists u k ∈ K such that

u k -ū 2 < 1 k and J(ū) + 1 2k u k -ū 2 2 > J(u k ). (2.8) Setting ρ k = u k -ū 2 and v k = (u k -ū)/ρ k , we can assume that v k ⇀ v in U 2 ;
if necessary, we select a subsequence. Then we prove that v ∈ C ū and J ′′ (ū)v 2 = 0. Because of (2.6), this is only possible if v = 0. With the help of (2.4) the contradiction is obtained from the identity v k 2 = 1; see [START_REF]Second order analysis for optimal control problems: Improving results expected from abstract theory[END_REF] for the details.

As a consequence of Theorem 2.2, we can not only prove that ū is the unique local minimum in a certain U 2 neighborhood. We are even able to show the nonexistence of other stationary points in such a neighborhood. Recall that ũ ∈ K is said to be a stationary point if

J ′ (ũ)(u -ũ) ≥ 0 for all u ∈ K.
(2.9)

Corollary 3. Under the assumptions of Theorem 2.2, there exists ε > 0 such that there is no stationary point ũ ∈ B 2 (ū; ε) ∩ K different from ū.

Assumption (2.6) has another consequence that was known up to now only in an U ∞ -neighborhood of ū. The result expresses some alternative formulation of second-order sufficient conditions that is useful for applications in the numerical analysis.

Theorem 2.4. Under the assumptions of Theorem 2.2, there exist numbers ε > 0, ν > 0 and τ > 0 such that

J ′′ (u)v 2 ≥ ν 2 v 2 2 ∀ v ∈ E τ ū and ∀ u ∈ K ∩ B 2 (ū; ε), (2.10
)

where E τ ū = {v ∈ cl 2 (S ū) : |J ′ (ū)v| ≤ τ v 2 }.

Application. An Elliptic Neumann Control Problem

In this section we study the optimal control problem

(P 1 ) min u∈K J(u),
where

J(u) = Ω L(x, y u (x)) dx + Γ l(x, y u (x), u(x)) dσ(x), (3.1)
K = {u ∈ L ∞ (Γ ) : α ≤ u(x) ≤ β for a.a. x ∈ Γ },
σ denotes the Lebesgue surface measure, -∞ < α < β < +∞, and y u is the solution of the following Neumann problem

-∆y + f (y) = 0 in Ω, ∂ ν y = u on Γ. (3.2)
We impose the following assumptions on the functions and parameters appearing in the control problem (P 1 ).

Assumption (N1):

Ω is an open, bounded and connected subset of R n , n ≥ 2, with Lipschitz boundary Γ and f : R

-→ R is a function of class C 2 such that f ′ (t) ≥ c o > 0 for all t ∈ R.
The reader is referred to [START_REF] Casas | Second order optimality conditions for semilinear elliptic control problems with finitely many state constraints[END_REF] for more general non-linear terms in the state equation.

Assumption (N2):

We assume that L :

Ω × R -→ R and l : Γ × R × R -→ R
are Carathéodory functions of class C 2 with respect to the second variable for L and with respect to the second and third variables for l, with L(•, 0) ∈ L 1 (Ω), l(•, 0, 0) ∈ L 1 (Γ ). For every M > 0 there exist functions ψ M ∈ L p(Ω), p > n/2, and φ M ∈ L q (Γ ), q > n -1, and a constant C M > 0 such that

                   ∂ j L ∂y j (x, y) ≤ ψ M (x), with j = 1, 2, ∂ j l ∂y j (x, y, u) ≤ φ M (x), with j = 1, 2, ∂ i+j l ∂u i ∂y j (x, y, u) ≤ C M , 1 ≤ i + j ≤ 2 and i ≥ 1
are satisfied for a.a. x ∈ Ω and every u, y ∈ R, with |y| ≤ M and |u| ≤ M . Moreover, for every ε > 0 there exists η > 0 such that for a.a. x ∈ Ω and all

u i , y i ∈ R, with i = 1, 2,      |y 2 -y 1 | ≤ η ⇒ ∂ 2 L ∂y 2 (x, y 2 ) - ∂ 2 L ∂y 2 (x, y 1 ) ≤ ε, |u 2 -u 1 | + |y 2 -y 1 | ≤ η ⇒ D 2 (y,u) l(x, y 2 , u 2 ) -D 2 (y,u) l(x, y 1 , u 1 ) ≤ ε.
Here D 2 (y,u) l(x, y, u) denotes the Hessian matrix of l with respect to the variables (y, u). We also assume the Legendre-Clebsch type condition

∃Λ > 0 such that ∂ 2 l ∂u 2 (x, y, u) ≥ Λ for a.a. x ∈ Γ and ∀y, u ∈ R. (3.3)
It is obvious that the usual quadratic integrands L(x, y) = 1 2 (y -y Ld (x)) 2 and l(x, y, u) = 1 2 (y -y ld (x)) 2 + Λ 2 u 2 satisfy Assumption (N2) if y Ld ∈ L p(Ω) and y ld ∈ L q (Γ ).

The hypothesis (3.3) is crucial for satisfying the assumptions (2.3) and (2.4). In §4 we will consider the case where (3.3) doe not hold.

On the state equation (2.1), the following result is known.

Theorem 3.1. Under the Assumption (N1), for every u ∈ L q (Γ ) the equation (3.2) has a unique solution y u ∈ H 1 (Ω) ∩ C( Ω). Furthermore, the mapping

G : L q (Γ ) -→ H 1 (Ω)∩C( Ω), defined by G(u) = y u , is of class C 2 . For elements u, v, v 1 and v 2 of L q (Γ ), the functions z v = G ′ (u)v and z v1v2 = G ′′ (u)(v 1 , v 2 )
are the solutions of the problems

Az + f ′ (y u )z = 0 in Ω, ∂ ν A z = v on Γ, (3.4)
and

Az + f ′ (y u )z + f ′′ (y u )z v1 z v2 = 0 in Ω, ∂ ν A z = 0 on Γ, (3.5)
respectively, where

z vi = G ′ (u)v i , i = 1, 2.
The proof of existence and uniqueness of a solution y u in H 1 (Ω) ∩ L ∞ (Ω) is standard; see, for instance, [START_REF] Casas | Boundary control of semilinear elliptic equations with pointwise state constraints[END_REF]. For the continuity of y u , the reader is referred to [START_REF] Haller-Dintelmann | Hölder continuity and optimal control for nonsmooth elliptic domains[END_REF] or [START_REF] Stampacchia | Problemi al contorno ellittici con dati discontinui dotati di soluzioni Hölderiane[END_REF]. As usual, the differentiability of G can be obtained from the implicit function theorem.

As a consequence of this theorem and the chain rule the next result follows.

Theorem 3.2. Assuming (N1) and (N2), then the mapping J : L ∞ (Γ ) -→ R, defined by (3.1), is of class C 2 . For all u, v, v 1 and v 2 of L ∞ (Γ ) we have

J ′ (u)v = Γ ϕ u + ∂l ∂u (x, y u , u) v dσ (3.6) J ′′ (u)(v 1 , v 2 ) = Ω ∂ 2 L ∂y 2 (x, y u ) -ϕ u f ′′ (y u ) z v1 z v2 dx + Γ ∂ 2 l ∂y 2 (x, y u , u)z v1 z v2 + ∂ 2 l ∂y∂u (x, y u , u)(v 1 z v2 + v 2 z v1 ) dσ + Γ ∂ 2 l ∂u 2 (x, y u , u)v 1 v 2 dσ, (3.7)
where

z vi = G ′ (u)v i , i = 1, 2, and ϕ u ∈ H 1 (Ω) ∩ C( Ω) is the solution of        -∆ϕ + f ′ (y u )ϕ = ∂L ∂y (x, y u ) in Ω, ∂ ν ϕ = ∂l ∂y (x, y u , u) on Γ. (3.8)
From the above expressions for J ′ (u) and J ′′ (u) and Assumption (N2) we deduce that J ′ (u) and J ′′ (u) can be extended to linear and bilinear forms, respectively, on L 2 (Γ ). Even more, there exist two constants M 1 > 0 and M 2 > 0 such that for every v, v 1 , v 2 ∈ L 2 (Γ ) and u ∈ K

|J ′ (u)v| ≤ M 1 v L 2 (Γ ) and |J ′′ (u)(v 1 , v 2 )| ≤ M 2 v 1 L 2 (Γ ) v 2 L 2 (Γ ) . (3.9)
This shows that (2.1) holds with U 2 = L 2 (Γ ) and U ∞ = L ∞ (Γ ). The most delicate issue in the proof of (2.2)-(2.4) is the verification of (2.3), which can be done with the help of the following lemma. Lemma 3.1 Let (X, Σ, µ) be a measure space with µ(X) < +∞. Suppose that

{g k } ∞ k=1 ⊂ L ∞ (X) and {v k } ∞ k=1 ⊂ L 2 (X) satisfy the assumptions -g k ≥ 0 a.e. in X, {g k } ∞ k=1 is bounded in L ∞ (X) and g k → g in L 1 (X). -v k ⇀ v in L 2 (X).
Then there holds the inequality

X g(x)v 2 (x) dµ(x) ≤ lim inf k→∞ X g k (x)v 2 k (x) dµ(x). (3.10)
The proof of this lemma can be obtained by an application of Egorov's theorem; see [START_REF]Second order analysis for optimal control problems: Improving results expected from abstract theory[END_REF]. To confirm (2.3) we apply Lemma 3.1 with X = Γ , µ = σ and

0 < Λ ≤ g k (x) = ∂ 2 l ∂u 2 (x, y u k (x), u k (x)) → g(x) = ∂ 2 l ∂u 2 (x, y u (x), u(x)) in L 1 (Γ ).
Finally, we apply Theorems 2.1 and 2.2 to the problem (P 1 ). Given ū ∈ K, we see that the cone of critical directions C ū defined in §2 can be expressed for the problem (P 1 ) in the form

C ū = {v ∈ L 2 (Γ ) : v(x) =    ≥ 0 if ū(x) = α ≤ 0 if ū(x) = β 0 if d(x) = 0 a.e. in Γ }, where d(x) = φ(x) + ∂l ∂u (x, ȳ(x), ū(x))
and ȳ = y ū and φ = ϕ ū denote the state and adjoint state associated to ū, respectively. It is not difficult to check that the regularity assumption stated in Theorem 2.1 is fulfilled by C ū. Then we have the following corollaries.

Corollary 3.1 Let the Assumption (N1) be satisfied and suppose that ū is a local minimum of (P 1 ) in the L ∞ (Γ ) sense. Then there holds

J ′ (ū)(u -ū) ≥ 0 for all u ∈ K and J ′′ (ū)v 2 ≥ 0 ∀v ∈ C ū. Conversely, if ū ∈ K obeys J ′ (ū)(u -ū) ≥ 0 ∀u ∈ K, (3.11) J ′′ (ū)v 2 > 0 ∀v ∈ C ū \ {0}, (3.12)
then there exist ε > 0 and δ > 0 such that

J(ū) + δ 2 u -ū 2 L 2 (Γ ) ≤ J(u) ∀u ∈ K ∩ B 2 (ū; ε). (3.13)
Let us underline that the mapping G is only differentiable in L q (Γ ) for q > n -1. Consequently, for all n ≥ 3, G is not differentiable in L 2 (Γ ). Moreover, the general nonlinear cost functional J is only differentiable in L ∞ (Γ ). Hence, for any dimension n, the classical theory of second order conditions would only assure the local optimality of ū in the L ∞ (Γ ) sense. In contrast to this, our result guarantees local optimality in the sense of L 2 (Γ ). Corollary 3.2 Under the assumption (N1) and (N2), there exists a ball B 2 (ū; ε) in L 2 (Γ ) such that there is no other stationary point in B 2 (ū; ε) ∩ K than ū. Moreover, there exist numbers ν > 0 and τ > 0 such that

J ′′ (u)v 2 ≥ ν 2 v 2 L 2 (Γ ) ∀ v ∈ C τ ū and ∀ u ∈ A ∩ B 2 (ū; ε), (3.14)
where A is a bounded open subset of L ∞ (Γ ) containing K and

C τ ū = {v ∈ L 2 (Γ ) : v(x) =    ≥ 0 if ū(x) = α ≤ 0 if ū(x) = β 0 if | d(x)| > τ
a.e. in Γ }.

In the above corollaries, B 2 (ū; ε) denotes the L 2 (Γ )-ball of radius ε centered at ū. Observe that the above cone C τ ū is not equal to the cone E τ ū defined in Theorem 2.4. However, if v ∈ C τ ū , then

|J ′ (ū)v| = Γ | d(x)v(x)| dx ≤ τ {x:| d(x)|≤τ } |v(x)| dx ≤ τ |Γ | v L 2 (Γ ) .
Thus, we have that

C τ ū ⊂ E τ Γ ū , with τ Γ = τ |Γ |.
Hence, Theorem 2.4 can be applied.

A Bang-Bang Control Problem

The reader is referred to [START_REF]Second order analysis for bang-bang control problems of pde[END_REF] for proofs and extensions of the results stated below.

Let Ω be an open and bounded domain in R n , n ≤ 3, with a Lipschitz boundary Γ . In this domain, we consider the following control problem

(P 2 )      min J(u) = Ω L(x, y u (x)) dx α ≤ u(x) ≤ β
where y u is the solution of the Dirichlet problem

-∆y + f (y) = u in Ω, y = 0 on Γ, (4.1)
-∞ < α < β < +∞ and L and f satisfy the following assumptions. Assumption (D1) The function f : R -→ R is of class C 2 and f ′ (t) ≥ 0 for every t ∈ R. Assumption (D2) The function L : Ω × R -→ R is measurable with respect to the first variable and of class C 2 with respect to the second. Moreover, L(•, 0) ∈ L 1 (Ω), and for all M > 0 there is a constant C L,M > 0 and a function ψ M ∈ L p(Ω) such that

∂L ∂y (x, y) ≤ ψ M (x), ∂ 2 L ∂y 2 (x, y) ≤ C L,M .
For every M > 0 and ε > 0 there exists δ > 0, depending on M and ε such that

∂ 2 L ∂y 2 (x, y 2 ) - ∂ 2 L ∂y 2 (x, y 1 ) < ε if |y 1 |, |y 2 | ≤ M, |y 2 -y 1 | ≤ δ, for a.a. x ∈ Ω.
Hereafter, we will denote

K = {u ∈ L ∞ (Ω) : α ≤ u(x) ≤ β for a.e. x ∈ Ω}.
For every u ∈ L p (Ω), with p > n/2, the state equation (4.1) has a unique solution y u ∈ H 1 0 (Ω) ∩ C( Ω). The proof of this result is a quite standard combination of Schauder's fixed point theorem and the L ∞ (Ω) estimates [START_REF] Stampacchia | Problemi al contorno ellittici con dati discontinui dotati di soluzioni Hölderiane[END_REF]. For the continuity of the solution in Ω see, for instance, [START_REF] Gilbarg | Elliptic Partial Differential Equations of Second Order[END_REF]Theorem 8.30]. Moreover, the mapping G :

L p (Ω) -→ H 1 0 (Ω) ∩ C( Ω), with G(u) = y u , is of class C 2 .
In the sequel, we will take p = 2 and we will denote by

z v = G ′ (u)v the solution of -∆z + f ′ (y u )z = v in Ω, z = 0 on Γ, (4.2)
where y u = G(u) is the state corresponding to u. As usual, we consider the adjoint state equation associated with a control u

   -∆ϕ + f ′ (y u )ϕ = ∂L ∂y (x, y u ) in Ω, ϕ = 0 on Γ, (4.3)
denoted by ϕ u . Because of the assumptions on L, we have that ϕ ∈ H 1 0 (Ω) ∩ C( Ω). Moreover, there exists M > 0 such that

y u ∞ + ϕ u ∞ ≤ M ∀u ∈ K. (4.4)
Under the above assumptions, the problem (P 2 ) has at least one solution ū with an associated state ȳ ∈ H 1 0 (Ω)∩C( Ω). The cost functional J : L 2 (Ω) -→ R is of class C 2 and the first and second derivatives are given by

J ′ (u)v = Ω ϕ u (x)v(x) dx, (4.5)
and

J ′′ (u)(v 1 , v 2 ) = Ω ∂ 2 L ∂y 2 (x, y u (x)) -ϕ u (x)f ′′ (y u (x)) z v1 (x)z v2 (x) dx, (4.6)
where

z vi = G ′ (v i ) are the solution of (4.2) for v = v i , i = 1, 2.
Any local solution ū satisfies the optimality system

-∆ȳ + f (ȳ) = ū in Ω, ȳ = 0 on Γ, (4.7)
   -∆ φ + f ′ (ȳ) φ = ∂L ∂y (x, ȳ) in Ω, φ = 0 on Γ, (4.8)
Ω φ(x)(u(x) -ū(x)) dx ≥ 0 ∀u ∈ K. (4.9)
From the last condition, we deduce as usual for a.a.

x ∈ Ω ū(x) = α if φ(x) > 0, = β if φ(x) < 0, and φ(x)    > 0 if ū(x) = α, < 0 if ū(x) = β, = 0 if α < ū(x) < β. (4.10)
The cone of critical directions associated with ū is defined by

C ū = {v ∈ L 2 (Ω) : v(x)    ≥ 0 if ū(x) = α ≤ 0 if ū(x) = β = 0 if φ(x) = 0 }
Then, the necessary second order condition satisfied is written in the form

J ′′ (ū)v 2 ≥ 0 ∀v ∈ C ū. (4.11)
For the above results the reader is referred to [START_REF] Casas | Second order optimality conditions for semilinear elliptic control problems with finitely many state constraints[END_REF] or [START_REF] Casas | Optimality conditions for a class of optimal control problems with quasilinear elliptic equations[END_REF], where similar cases were studied. Let us remark that in the case where the set of zeros of φ has a zero Lebesgue measure, then ū(x) is either α or β for almost all points x ∈ Ω, i.e. ū is a bang-bang control. Moreover, in this case, C ū = {0}, therefore (4.11) does not provide any information. Consequently, it is unlikely that the sufficient second order conditions could be based on the set C ū. To overcome this drawback we are going to increase the set C ū. For every τ ≥ 0 we define

C τ ū = {v ∈ L 2 (Ω) : v(x)    ≥ 0 if ū(x) = α ≤ 0 if ū(x) = β = 0 if | φ(x)| > τ

}

It is obvious that C 0 ū = C ū. An example due to Dunn [START_REF] Dunn | Second-order optimality conditions in sets of L ∞ functions with range in a polyhedron[END_REF] proves that, in general, the second order condition based on the cone C ū is not sufficient for the local optimality. Before analyzing (P 2 ), let us take a look on its Tikhonov regularization. Then, we have

J ′ Λ (u)v = Ω (ϕ u + Λu)v dx and J ′′ Λ (u)(v 1 , v 2 ) = Ω ∂ 2 L ∂y 2 (x, y u ) -ϕ u ∂ 2 f ∂y 2 (x, y u ) z v1 z v2 dx + Λ Ω v 1 v 2 dx.
Now, we apply Theorem 2.2 to (P 2,Λ ) and we obtain the following result.

  For any Λ > 0, let us consider the problem(P 2,Λ ) min u∈K J Λ (u) = Ω L(x, y u (x)) dx + Λ 2 Ω u 2 (x) dx.

Theorem 4 . 1 .

 41 Let ū ∈ K satisfy that J ′ Λ (ū)(u -ū) ≥ 0 ∀u ∈ K and J ′′ Λ (ū)v 2 > 0 ∀v ∈ C ū \ {0}.

⋆ This work was partially supported by the Spanish Ministerio de Economía y Competitividad under project MTM2011-22711.

Then, there exists δ > 0 and ε > 0 such that

In the above theorem and hereafter, B 2 (ū; ε) denotes the L 2 (Ω)-ball of center at ū and radius ε. Now, invoking Theorem 2.4 and observing that C τ ū ⊂ E τ Ω ū for τ Ω = |Ω|τ , we get the following theorem.

Then, the following assumptions are equivalent

where

In the case Λ = 0, Dunn's example shows that 1 is not enough, in general, to assure the local optimality of ū. We will see below that 2 does not hold for Λ = 0. Then, it remains to analyze if the assumption 3 is enough for the local optimality of ū when Λ = 0. The next theorem proves that it is sufficient. Theorem 4.3. Let us assume that ū is a feasible control for problem (P 2 ) satisfying the first order optimality conditions (4.7)-(4.9) and suppose that there exist δ > 0 and τ > 0 such that

where z v = G ′ (ū)v is the solution of (4.2) for y = ȳ. Then, there exists ε > 0 such that

with z u-ū = G ′ (ū)(u -ū).

Corollary 4.1 Under the hypotheses of Theorem 4.3, there exists ε > 0 such that

We finish by showing that the statement 2 of Theorem 4.2 does not hold for Λ = 0. Indeed, let us assume that it holds. Then, a simple modification of the proof of Theorem 4.3, see [START_REF]Second order analysis for bang-bang control problems of pde[END_REF], leads to the inequality

for some ν > 0 and ε > 0. Then, ū is a solution of the problem

The Hamiltonian of this control problem is given by

From the Pontryagin's principle we deduce

H(x, ȳ(x), t, φ(x)) for almost all x ∈ Ω.

However, invoking (4.10) we obtain that this is a contradiction to the following facts that can be easily checked


Michel C Delfour
email: delfour@crm.umontreal.ca

André Garon
email: andre.garon@polymtl.ca

Quadratic ODE and PDE models of drug release kinetics from biodegradable polymers

Keywords: Drug release models, biodegradable polymers, paclitaxel

In order to achieve prescribed drug release kinetics over long therapeutic periods, bi-phasic and possibly multi-phasic releases from blends of biodegradable polymers are currently envisioned. The modelling of drug release in the presence of degradation of the polymer matrix and surface erosion is quite complex. Yet, simple reliable mathematical models validated against experimental data are now available to help in classifying neat polymers and in predicting the release dynamics from polymer blends. In this paper, we survey a two-parameter quadratic ODE model that has been validated against experimental data for the release of paclitaxel from a broad range of biodegradable polymers and a quadratic semi-permeable membrane PDE model that mimics the ODE model and could readily be extended to drug eluding stents.

Introduction

Stents are used in interventional cardiology to keep a diseased vessel open after angioplasty. This procedure is known to damage the endothelium at the insertion site and thus to favour the occurrence of in-stent restenosis through the proliferation of smooth muscle cells (SMC) within the vessel lumen. To control the abnormal behaviour of SMC, stents are coated with polymers that slowly release drug through diffusion into the vessel wall (drug-eluting stents or DES). These drugs are designed to control the rate of mitosis of SMC until the regeneration of the endothelium. The reader is referred to T. Kataoka et als [START_REF] Kataoka | 7-Hexanoyltaxol-Eluting Stent for Prevention of Neointimal Growth: An Intravascular Ultrasound Analysis From the Study to Compare REstenosis rate between QueST and QuaDS-QP2 (SCORE)[END_REF] in 2002 and Joner et als [START_REF] Joner | Pathology of drug-eluting stents in humans -delayed healing and late thrombotic risk[END_REF] in 2006 for a fairly well-documented account of DES for the prevention of neointimal growth (see, for instance, [15, Figure 1,p. 1791]).

If endothelial cells do not recover to effectively control the proliferation of SMC's, a sustained dose will be required over the therapeutic period and even forever. In order to achieve prescribed drug release kinetics the current design ⋆ Centre de recherches mathématiques et Département de mathématiques et de statistique, Université de Montréal, CP 6128, succ Centre-ville, Montréal (Qc), Canada H3C 3J7. ⋆⋆ Département de Génie mécanique, École Polytechnique de Montréal, C.P. 6079, succ. Centre-ville, Montréal (Qc), Canada H3C 3A7.

strategies focus on bi-phasic and possibly multi-phasic releases from blends of biodegrable polymers (see, for instance, Batycky et al [START_REF] Batycky | A theoretical model of erosion and macromolecular drug release from biodegrading microspheres[END_REF] in 1997) to achieve specific drug release kinetics profiles over long therapeutic windows.

Recently, Lao and Venkatraman [START_REF] Lao | Adjustable paclitaxel release kinetics and its efficacy to inhibit smooth muscle cells proliferation[END_REF] and Lao, Venkatraman, and Peppas [START_REF] Lao | Modeling of drug release from biodegradable polymer blends[END_REF] have proposed a semi-empirical model to predict the release profile of paclitaxel from three neat polymer matrices: PCL (Polycaprolactone), PLGA (dl-lactideco-glycolide) and PLGAPEG (PLGA with polyethylene glycol). They are representative of a broad family of biodegradable polymers ranging from hydrophobic to hydrophilic. In hydrophilic polymers the internal bounds between the chains are weakened and this adds to the surface erosion phenomenon. The drug release mechanism within a polymer matrix depends on many factors such as the affinity of the drug with the surrounding medium (water). Specifically, paclitaxel is hydrophobic and this might explain the fact that some of the drug blended into the polymer matrix is not released and cannot participate to the treatment of the disease wall. This is a difficult subject. The main criticism expressed in [START_REF] Lao | Modeling of drug release from biodegradable polymer blends[END_REF] of available models for drug release from eroding surfaces is that they fail to faithfully reproduce experimental data for highly degradable polymers (the S-curve behaviour). The reader is referred to the introduction of the paper of Lao et als [START_REF] Lao | Modeling of drug release from biodegradable polymer blends[END_REF] for a comprehensive review of the literature.

A quick look at the paclitaxel release profiles suggests two types of release: S-curve type and exponential type. S-curve behaviours are similar to the ones encountered in the study of the logistic equation of populations. In [START_REF] Blanchet | Quadratic models to fit experimental data of Paclitaxel release kinetics from biodegradable polymers[END_REF] we introduced a simple two-parameter Ordinary Differential Equation (ODE) model that completely describes the paclitaxel release profiles from neat PCL, PLGA PEG, and PLGA polymer matrices. This model describes with greater accuracy the drug-release than the semi-empirical model of Lao et als [START_REF] Lao | Modeling of drug release from biodegradable polymer blends[END_REF] using 5 to 8 parameters.

The simplicity of our model for such a broad range of polymers indicates that somehow the quadratic structure captures the complex microphysics and chemistry of the release and degradation processes. Using a purely mathematical intuition to modelling, we have introduced in [6] a time-space three dimensional partial differential equation (PDE) model of the paclitaxel release that mimics the ODE model. The film of neat polymer is modelled as a thin flat domain whose polymer/medium interface is a quadratic semi-permeable membrane with a concentration jump at the interface.

In this approach, the diffusion process through a semi-permeable membrane is modelled as a diffusion through an interface with cracks (not to be confused with holes) where the rate of transfer of the product is proportional to the size of the concentration jump across the interface. Since the cracks have zero surface, their size is measured in terms of the mathematical notion of capacity. What is very nice about this approach is that it is based on a mathematically well documented linear model coming from the study of the Neumann sieve by Damlamian [START_REF] Damlamian | Le problème de la passoire de Neumann, (French) [The Neumann sieve problem[END_REF] in 1985. It provides a variational formulation and a mathematically tractable approach to the asymptotic analysis of a punctured membrane as the size of the holes goes to zero while preserving a strictly positive capacity that accounts for the diffusion of the drug through the cracks 1 Adding the non-linearity captures the effect of the internal degradation of the polymer by making the rate of mass transfer proportional to the size of the concentration jump across the interface.

Our approach is different from others in the literature since it deals with the nonlinearity through a quadratic condition at the interface between the polymer and the surrounding medium instead of using a time-dependent or a nonlinear diffusion. This model can be seen as a first step towards a three dimensional modelling of the release of paclitaxel from drug eluding stents coated with biodegradable polymers. It is capable of covering a wide range of biodegradable polymers potentially including the ones for which an incomplete release is experimentally observed (recall that the paclitaxel is hydrophobic).

To complete the experimental approach to this modelling, the next step would be to set up an experimental benchmark to check if the model and the mathematical assumptions on the coefficients of the model are realistic. The validation of such a model would improve the modelling of the drug release part of the global three-dimensional model of a blood vessel incorporating the lumen, the blood, the aggregated wall, and the coated stent (cf. for instance, [START_REF] Delfour | Modeling and design of stents to optimize the effect of the dose[END_REF]) and the subsequent studies of the effect of the pattern of the stent in [START_REF] Bourgeois | General patterns and asymptotic dose in the design of coated stents[END_REF] and the effect of the pulsative nature of the blood in [START_REF] Delfour | New equations for the dose under pulsative/periodic conditions in the design of coated stents[END_REF]. Such global studies are important to determine the set of features in the modelling of the blood vessel and of the stent that should be retained in the design of the stent and the drug release dynamics.

ODE model and gradient flow interpretation

In the previous paper [START_REF] Blanchet | Quadratic models to fit experimental data of Paclitaxel release kinetics from biodegradable polymers[END_REF] we have shown an excellent fit between experimental release data [START_REF] Lao | Adjustable paclitaxel release kinetics and its efficacy to inhibit smooth muscle cells proliferation[END_REF] of paclitaxel from biodegradable neat polymers and a twoparameter quadratic ODE model of the Riccati type. We briefly recall this model.

Given an initial mass M 0 > 0 of drug uniformly impregnated into a polymeric matrix, denote by M (t) > 0 the mass of drug released outside the polymer as a function of the time t > 0. Denote by M ∞ , 0 ≤ M ∞ ≤ M 0 , the asymptotic mass of the drug released. The ODE model was chosen of the form

dM dt (t) = h(M (t)), t > 0, M (0) = 0, (2.1)
for some quadratic right-hand side

h(M ) def = A 1 (M ∞ -M (t)) + A 2 (M ∞ -M (t)) 2 (2.2) such that M ′ (0) = (A 1 + A 2 M ∞ ) M ∞ > 0. By introducing the normalized released mass m(t) def = M (t)/M 0 , (2.3)
we get the following quadratic ODE model

dm dt (t) = A 1 + A 2 M 0 M ∞ M 0 -m(t) M ∞ M 0 -m(t) , m(0) = 0. (2.4)
Assuming that the ratio 0 < M ∞ /M 0 ≤ 1 is known, the model is completely specified by the two parameters A 1 and A 2 M ∞ . When A 2 = 0, the model is linear; when A 2 = 0, the right-hand side is of the form

h(m) def = A 2 M 0 (m 2 -m)(m 1 -m), m 1 def = M ∞ M 0 , m 2 def = A 1 + A 2 M ∞ A 2 M 0 .
It was shown in [START_REF] Blanchet | Quadratic models to fit experimental data of Paclitaxel release kinetics from biodegradable polymers[END_REF] that the following four cases can occur under the conditions m(0) = 0 and m ′ (0

) = A 1 + A 2 M ∞ > 0: Case 1) (True S type) A 1 > 0, A 2 < 0, and -m 1 < 1 2 A 1 A 2 M 0 (that is, -m 1 < m 2 < 0), (2.5)
with solution

m(t) = m 1 m 2 1 -e -A1 t m 2 -m 1 e -A1 t
for which the point of inflexion occurs at time

t c = -(log(-m 2 /m 1 ))/A 1 > 0; Case 2) (S type) A 1 > 0, A 2 < 0, and 1 2 A 1 A 2 M 0 ≤ -m 1 (that is, m 2 ≤ -m 1 ), (2.6)
with the solution and the point of inflexion

m(t) = m 1 m 2 1 -e -A1 t m 2 -m 1 e -A1 t , t c = -(log(-m 2 /m 1 ))/A 1 ≤ 0; Case 3) (Exponential type) A 1 ≥ 0 and A 2 > 0 (that is, m 2 ≥ 1), (2.7)
with the solution and the blow up time

       m(t) = m 1 m 2 1 -e -A1 t m 2 -m 1 e -A1 t , t c = - log(m 2 /m 1 ) A 1 < 0, for A 1 > 0 since m 2 > m 1 ,        m(t) = m 1 A 2 M ∞ t 1 + A 2 M ∞ t , t c = - 1 A 2 M ∞ < 0, for A 1 = 0 since m 2 = m 1 ;
(2.8)

Case 4) (True exponential) A 1 > 0 and A 2 = 0 with the solution

m(t) = M ∞ M 0 1 -e -A1 t , t c = -∞.
The generic behaviours of the solution m(t) in the above cases are illustrated in Figures 1 and2 and the parameters tabulated in Table 1 of [START_REF] Delfour | Drug release kinetics from biodegradable polymers via partial differential equations models[END_REF].

The ODE model has an interesting gradient flow interpretation by introducing the function

E(m) def = 1 2 A 1 M ∞ M 0 -m 2 + 1 3 A 2 M 0 M ∞ M 0 -m 3 (2.9)
with gradient (derivative)

E ′ (m) = -A 1 M ∞ M 0 -m -A 2 M 0 M ∞ M 0 -m 2 (2.10)
and Hessian (second order derivative)

E ′′ (m) = A 1 + 2 A 2 M 0 M ∞ M 0 -m . (2.11)
The ODE can now be rewritten in the form of a gradient flow equation

dm dt (t) + E ′ (m(t)) = 0, m(0) = 0. (2.12)
This is the continuous version of a steepest descent method to minimize the functional E. So, it is expected that starting from m(0) = 0 with m ′ (0) > 0 the asymptotic value m 1 of the solution of the ODE (2.12) would achieve a local minimum of E(m). To do that, we compute the second derivative of E under the assumption that A 1 ≥ 0 and m ′ (0) > 0 which is equivalent to

E ′ (0) = -(M ∞ /M 0 )[A 1 + A 2 M ∞ ] < 0.
It turns out that in all cases except the second part of case 3), m 1 is a local minimum of E(m). The exception corresponds to a point of inflection that can be changed into a global minimum by modifying the function

E to E(m) = (A 2 M 0 /3) |M ∞ /M 0 -m| 3 .
3 PDE model of quadratic semi-permeable membranes

Equations in the polymer and the surrounding medium

The experimental benchmark of [START_REF] Lao | Adjustable paclitaxel release kinetics and its efficacy to inhibit smooth muscle cells proliferation[END_REF] is contained in a vial. The polymer film is deposited flat at the bottom of the vial and the vial is filled with a fluid that we shall call the surrounding medium (see Figure 1).

∂c m ∂t = div (D m ∇c m ) in Ω m (3.2)
with constant diffusion constants D p and D m and initial conditions

c p (x, 0) = c 0 (x) = M 0 /|Ω p | in Ω p , c m (x, 0) = 0 in Ω m , (3.3)
where |Ω p | is the volume of Ω p . Assume that the experimental set up is closed:

D p ∂c p ∂n p = 0eq.constriantΓ p \Γ int D m ∂c m ∂n m = 0 on Γ m \Γ int , (3.4)
where the unit normals n p and n m are exterior to the respective domains Ω p and Ω m . Assume that there is no loss of product: this yields the (affine) constraint

∀t ≥ 0, M 0 def = Ωp c 0 (x) dx = Ωp c p (x, t) dx + Ωm c m (x, t) dx, (3.5)
where M 0 is the total mass of product. By integrating (3.1) over Ω p and (3.2) over Ω m and by using the constraint (3.5), we get It remains to specify the conditions at the interface Γ int .

⇒ 0 = Ωp ∂c p ∂t (x, t) dx + Ωm ∂c m ∂t (x, t) dx = Ωp div (D p ∇c p )(x, t) dx + Ωm div (D m ∇c m )(x,

Conditions at the interface

In order to incorporate the microphysics taking place in the thin film of polymer, it is assumed that the interface behaves as a semi-permeable membrane with micro fissures through which the drug diffuses into the surrounding medium. Many empirical and theoretical models of such membranes have been studied in the literature and in different contexts. One mathematically interesting model of a semi-permeable membrane is to assume that the interface is a membrane punctured with small holes whose size goes to zero while preserving a strictly positive capacity2 in the limiting process. In other words the membrane is fissured or cracked and the drug diffuses through the cracks. This problem has been studied from the mathematical point of view under the name of the Neumann sieve by A. Damlamian [START_REF] Damlamian | Le problème de la passoire de Neumann, (French) [The Neumann sieve problem[END_REF] in 1985. From the physical point of view, it can be assimilated with a semi-permeable membrane.

In this section we consider an evolution equation of the form

∂c ∂t (t) + A(c(t)) = 0, c(0) = M 0 /|Ω p | χ Ωp , (3.8)
where the operator A is now quadratic in c(t). Since the domain Ω p is thin, it is reasonable to put the nonlinearity at the interface Γ int rather than on Ω p via a diffusion coefficient D p (c) that depends on c:

- d dt Ωp c p (t) dx = d dt Ωm c m (t) dx = Γint k 1 + k 2 |Ω p | M 0 |c p (t) -c m (t)| (c p (t) -c m (t)) dΓ (3.9)
for some constant k 2 . Note that we have introduced a scaling by the initial concentration of product M 0 /|Ω p | of the drug so that k 1 and k 2 are parameters of the same physical dimension. Now consider the (cubic) functional

E(v) def = 1 2 Ωp D p |∇v p | 2 dx + 1 2 Ωm D m |∇v m | 2 dx + Γint 1 2 k 1 |v p -v m | 2 + 1 3 k 2 |Ω p | M 0 |v p -v m | 3 dΓ (3.10) v p def = v| Ωp , v m def = v| Ωm (3.11) defined on the space H 1 (Ω p ∪ Ω m ) with a crack Γ int = Γ p ∩ Γ m in Ω p ∪ Ω m along which the function v can have a jump discontinuity [v] = v m -v p .
This convex non quadratic variational formulation is similar to the T 4 radiation law for the temperature T of a radiating body in free space (cf., for instance, [START_REF] Delfour | Approximation of nonlinear problems associated with radiating bodies in space[END_REF]). We do not impose the continuity of the concentrations at the interface. Taking into account the constraint on the total mass of product, we look for a solution c(t) at time t > 0 in the affine subspace

V pm M0 def =                c ∈ H 1 (Ω p ∪ Ω m ) : Ωp∪Ωm c(x) dx = M 0 , if k 2 = 0, k 1 > 0,      c ∈ H 1 (Ω p ∪ Ω m ) : Ωp∪Ωm c(x) dx = M 0 v p -v m ∈ L 3 (Γ int )      , if k 2 > 0, of H 1 (Ω p ∪ Ω m ). In the first case Ωp |∇v p | 2 dx + Ωm |∇v m | 2 dx + Γint |v p -v m | 2 dΓ 1/2 (3.12)
is an equivalent norm on V pm M0 ; in the second case

Ωp |∇v p | 2 dx + Ωm |∇v m | 2 dx 1/2 + Γint |v p -v m | 3 dΓ 1/3 (3.13)
is an equivalent norm on V pm M0 (cf., for instance, [START_REF] Delfour | Approximation of nonlinear problems associated with radiating bodies in space[END_REF] with the T 4 radiation law for the temperature T in free space).

The directional derivative of E is

dE(u; v) = Ωp D p ∇u • ∇v dx + Ωm D m ∇u • ∇v dx + Γint k 2 |Ω p | M 0 |u p -u m | (u p -u m ) (v p -v m ) + k 1 (u p -u m ) (v p -v m ) dΓ (3.14) u p def = u| Ωp , u m def = u| Ωm , v p def = v| Ωp , v m def = v| Ωm . (3.15)
We are interested in the stationary points c = (c p , c m ) ∈ V M0 of E that are the solutions of the variational equation

∃c ∈ V pm M0 , dE(c; v) = 0, ∀v ∈ V pm 0 , (3.16
)

V pm 0 def =                c ∈ H 1 (Ω p ∪ Ω m ) : Ωp∪Ωm c(x) dx = 0 , if k 2 = 0, k 1 > 0,      c ∈ H 1 (Ω p ∪ Ω m ) : Ωp∪Ωm c(x) dx = 0 v p -v m ∈ L 3 (Γ int )      , if k 2 > 0.
(3.17)

Again dE(c; v) = 0 for all constant functions v and V pm 0 can be replaced by

H 1 (Ω p ∪ Ω m ): ∃c ∈ V pm M0 , dE(c; v) = 0, ∀v ∈ H 1 (Ω p ∪ Ω m ).
It yields a complete set of conditions at the interface and the following system of equations div (D p ∇c p ) = 0 in Ω p , div (D p ∇c m ) = 0 in Ω m (3.18) ) on the temperature of a radiating body (cf., for instance, [START_REF] Delfour | Approximation of nonlinear problems associated with radiating bodies in space[END_REF]). The thin layer of polymer behaves as a nonlinear semi-permeable membrane. The second order directional derivative of E is

D p ∂c p ∂n p + k 2 |Ω p | M 0 |c p -c m | (c p -c m ) + k 1 (c p -c m ) = 0 on Γ int (3.19) D m ∂c m ∂n m -k 2 |Ω p | M 0 |c p -c m | (c p -c m ) + k 1 (c p -c m ) = 0 on Γ int (3.
d 2 E(u; v; w) = Ωp D p ∇w • ∇v dx + Ωm D m ∇w • ∇v dx + Γint 2 k 2 |Ω p | M 0 |u p -u m | + k 1 (w p -w m ) (v p -v m ) dΓ (3.23) ⇒ d 2 E(u; v; v) = Ωp D p |∇v| 2 dx + Ωm D m |∇v| 2 dx + Γint 2 k 2 |Ω p | M 0 |u p -u m | + k 1 |v p -v m | 2 dΓ. (3.24)
Since E is a cubic functional, local minima and local maxima can both occur depending on the signs and magnitudes of the constants k 1 and k

2 . A local minimum u ∈ V pm M0 is characterized by ∀v ∈ V pm 0 dE(u; v) = 0 and ∀0 = v ∈ V pm 0 d 2 E(u; v; v) > 0 and a local maximum u ∈ V pm M0 by ∀v ∈ V pm 0 dE(u; v) = 0 and ∀0 = v ∈ V pm 0 d 2 E(u; v; v) < 0.
Going back to the evolution equation (3.8) using the above conditions at the interface, we get the following system of equations

∂c p ∂t = div (D p ∇c p ) in Ω p , ∂c m ∂t = div (D m ∇c m ) in Ω m c p (x, 0) = M 0 /|Ω p | χ Ωp (x) in Ω p , c m (x, 0) = 0 in Ω m D p ∂c p ∂n p + k 2 |Ω p | M 0 |c p -c m | (c p -c m ) + k 1 (c p -c m ) = 0 on Γ int D m ∂c m ∂n m -k 2 |Ω p | M 0 |c p -c m | (c p -c m ) + k 1 (c p -c m ) = 0 on Γ int D p ∂c p ∂n p = 0 on Γ p \Γ int , D m ∂c m ∂n m = 0 on Γ m \Γ int Ωp c p dx + Ωm c m dx = M 0 .
(3.25)

The nonlinear condition on

Γ int D m ∂c m ∂n m = k 2 |Ω p | M 0 |c p -c m | (c p -c m ) + k 1 (c p -c m ) = k 2 |Ω p | M 0 |c p -c m | + k 1 k(c) (c p -c m )
says that k(c) is an affine function of the size of the jump. This means that the rate of transfer of the product across the interface is large when the absolute value of the concentration jump is large. Assuming that k 1 ≥ 0, when k 2 > 0 it decreases to k 1 when the size of the jump goes to zero; when k 2 < 0 it increases to k 1 when the size of the jump goes to zero.

Remark 1. When k 2 > 0, it would not be appropriate to remove the absolute value on c p -c m in the term k ′ of the previous identity. This would give the expression

D m ∂c m ∂n m = k 2 |Ω p | M 0 (c p -c m ) + k 1 k ′ (c) (c p -c m ),
where, if the size of the jump is large, k ′ (c) > 0 is large, ∂c m /∂n m > 0 is large, and the diffusion of product would be from the medium to the polymer even when c p > c m , that is, when the concentration in the polymer is larger than the one in the medium. However, it is interesting to note that various behaviours can be modelled by replacing

Fig. 1 .

 1 Fig. 1. The polymer film and the surrounding fluid in the vial.

  20) D p ∂c p ∂n p = 0 on Γ p \Γ int , D m ∂c m ∂n m = 0 on Γ m \Γ int (3.21) Ωp c p dx + Ωm c m dx = M 0 . (3.22) From the mathematical viewpoint, the condition involving |c p -c m | (c p -c m ) is the analogue of the condition |T -T m | 3 (T -T m ) (usually written (T -T m ) 4

  |c p -c m | by the plus [c p -c m ] + = max{0, c p -c m } or the minus [c p -c m ] -= max{0, -(c p -c m )} functions or introducing a threshold θ > 0 max{|c p -c m | -τ, θ}.

  h = |Ω p |/|Γ int | is the thickness of the polymer. This would correspond to A 1 = k 1 /h and A 2 = k 2 /h in the ODE model. The thickness h is an important parameter : the thinner the polymer the faster the release. If k 1 and k 2 are constants, m m can be normalized through the change of variable t → τ = t/h.

  3.3 Relation between the PDE and the ODE modelsSince |Ω p | is much smaller than |Ω m |, this last equation is related to the quadratic ODE model by making the same assumptions on the concentrations on Γ int as in the previous section:c p (x, t) ≃ 1 |Ω p | Ωp c p (x, t) dx and c m (x, t) ≃ 1 |Ω m | Ωm c m (x, t) dx (3.26) ⇒ c p (x, t) -c m (x, t) ≃ 1 |Ω p | [M 0 -M m (t)] ,

	where
								M m (t)	def =	c m (x, t) dx	(3.27)
								Ωm
	is the mass released at time t in the medium and
		dM m dt	(t) =	|Γ int | |Ω p |	k 1 +	k 2 M 0	|M 0 -M m (t)| (M 0 -M m (t))	(3.28)
	⇒	dm m dt	(t) =	1 h	[k

1 + k 2 |1 -m m (t)|] (1 -m m (t)), m m (t)

See also the more recent comprehensive paper [4, Theorem 5.5] using the very nice theory of periodic unfolding.

The capacity of a set is a mathematical notion. For instance a finite segment in the plane has zero area but finite capacity. Roughly speaking, the capacity is a "measure" of the cracks.