Xabier Iturbe
email: xabier.iturbe@arm.com

Didier Keymeulen
email: didier.keymeulen@jpl.nasa.gov

Patrick Yiu
email: pyiu@mit.edu

Daniel Berisford

Robert Carlson

Kevin Hand

Emre Ozer
email: emre.ozer@arm.com

⋆ Xabier Iturbe

On the Use of System-on-Chip Technology in Next-Generation Instruments Avionics for Space Exploration

Keywords: Fault-Tolerance, Avionics, System-on-Chip Integration, ARM Processor, Signal Processing

technology enables integrating all the functionality required to control and process science data delivered by space instruments in a single silicon chip (e.g., microprocessor + programmable logic). This chapter discusses the implications of using this technology in deep-space exploration avionics, namely in the next generation of NASA science instruments that will be used to explore our Solar system. We present here our experience at the NASA Jet Propulsion Laboratory (JPL) using Xilinx Zynq SoC devices to implement the data processing of a Fourier transform spectrometer, namely the Compositional InfraRed Imaging Spectrometer (CIRIS). Besides, we also discuss the different fault-tolerance techniques that have been implemented in the CIRIS controller SoC to deal with harsh radiation conditions prevailing in deep-space environments.

Introduction

Hybrid System-on-Chip (SoC) devices that embed the most energy efficient processor (ARM cores [START_REF] Blem | A Detailed Analysis of Contemporary ARM and x86 Architectures[END_REF]) and the latest and most powerful FPGA architecture (Xilinx 7-series [START_REF] Mehta | Xilinx 7 Series FPGAs: The Logical Advantage[END_REF]) into a single chip (Xilinx Zynq [START_REF] Inc | Zynq-7000 All Programmable SoC: Technical Reference Manual[END_REF]) promise new opportunities due to the performance, power consumption, weight and volume benefits they bring. This is especially relevant for building more capable space avionics. Currently most of these systems combine programmable logic and processors as separate components distributed along one or several PCB board(s), which results in power consumption overheads and larger volume to be put into space [START_REF] Wiens | The ChemCam Instrument Suite on the Mars Science Laboratory (MSL) Rover: Body Unit and Combined System Tests[END_REF][START_REF] Petrick | Adapting the SpaceCube v2.0 Data Processing System for Mission-Unique Application Requirements[END_REF]. Besides, currently existing space-grade processors (e.g., RAD750 [START_REF]RAD750 Radiation-Hardened PowerPC Microprocessor[END_REF]) are not suitable to be used in the next-generation spacecraft computing platforms because they do not provide sufficient performance and energy efficiency [START_REF] Doyle | High Performance Spaceflight Computing (HPSC) Next-Generation Space Processor (NGSP): A Joint Investment of NASA and AFRL[END_REF]. As a result, NASA and other space agencies have approached ARM and SoC technology, hoping to pave the way for future space exploration missions that are becoming ever more performance demanding.

Despite the fact that currently there are no space-qualified SoC parts, NASA is testing commercial Xilinx Zynq SoC devices in the International Space Station (ISS) as well as in precursor CubeSats operating in Low Earth Orbit (LEO), where the exposure to radiation is limited.

In view of a potential radiation-hardened SoC device that might be ready to fly in deep-space missions in the near to mid future, JPL and ARM have partnered together to develop a SoC platform to be used as a research vehicle for powering next-generation flight instruments intended to be used in NASA deep-space missions. Presently this platform, called APEX-SoC (APEX stands for Advanced Processor core for space EXploration), is being prototyped using a commercial Xilinx Zynq device. The APEX-SoC includes a generic and adaptable infrastructure that provides support for hardware and software based science processing. More specifically, the data acquisition and processing proper to each science instrument is to be implemented as a collection of "custom software and hardware applications" that are encapsulated by the APEX-SoC infrastructure and run on the Zynq's on-chip ARM processor and reside on the Zynq's FPGA fabric. Besides the infrastructure itself, the APEX-SoC includes a set of Radiation Hardened By Design (RHBD) features to protect the instrument-dependent modules implemented on the FPGA fabric from harsh space radiation. In connection with this, we are currently carrying out two research efforts to create a space-grade ARM processor that could potentially replace commercial ARM processors embedded in future radiation-tolerant SoC devices. First, we are conducting a thorough soft-error analysis of the ARM Cortex-R5 microprocessor, which is currently used in terrestrial safety-critical real-time applications, to identify the most vulnerable parts in the micro-architecture of this processor, analyze what level of protection is required for these vulnerable parts (e.g., detection only, correction only or hybrid), and then decide how to achieve this level of protection. Secondly, we are designing a Cortex-R5 based fail-operational Triple Core Lock-Step ARM processor (TCLS-ARM) with the capability to recover from errors within microseconds [START_REF] Iturbe | A Triple Core Lock-Step (TCLS) ARM Cortex-R5 Microprocessor for Safety-Critical and Ultra-Reliable Applications[END_REF].

This chapter describes the first prototype of the APEX-SoC platform implemented on the Zynq SoC and presents an illustrative case-study drawn from the JPL Compositional Infrared Imaging Spectrometer (CIRIS) [START_REF] Berisford | Thermal Testing of the Compositional Infrared Imaging Spectrometer (CIRIS)[END_REF], which has been proposed to be used in icy moons, such as Jupiter's moon Europa [START_REF] Carlson | The Compositional InfraRed Interferometric Spectrometer (CIRIS) for Assessing the Habitability of Europa[END_REF]. The remainder of this chapter is as follows. Section 2 introduces the SoC technology and its use in space missions so far. Section 3 describes the APEX-SoC platform, and then section 4 presents a case-study where the JPL CIRIS spectrometer data processing is implemented on this platform. Section 5 summarizes the implementation, performance and irradiation results that have been collected so far and, finally, section 6 concludes the chapter and points out to future work.

System-on-Chip Technology and its Use in Space Exploration Avionics

Despite miniaturized SoC technology is very convenient for space, where every gram of mass launched involves enormous costs, it is currently designed for and used by consumer terrestrial applications, where a single device with very low power consumption has found a niche in the network and telecommunication markets. Commercial SoCs have developed very advanced computation capabilities in consumer electronics that is continuously demanding more powerful devices and applications. One example of the sophistication degree achieved by commercial SoCs is the Xilinx Zynq-UltraScale+ MPSoC that is scheduled for release in early 2016 [START_REF] Inc | UltraScale Architecture and Product Overview[END_REF]. This will include an ARM Cortex-A53 high-performance 64-bit processor, an ARM Cortex-R5 real-time processor and a Xilinx UltraScale FPGA architecture. Current SoC devices available in the market typically include at least one processor and an FPGA fabric. Since ARM cores are the standard processors used in all SoCs, the difference between them comes from the FPGA fabric they use. This fabric embeds routing resources, programmable logic, DSP and RAM blocks together with the memory cells to store their configuration.

Although the current use of SoCs is largely limited to terrestrial applications, space agencies consider this technology could be an alternative to overcome the current performance crisis seen in the space sector [START_REF] Doyle | High Performance Spaceflight Computing (HPSC) Next-Generation Space Processor (NGSP): A Joint Investment of NASA and AFRL[END_REF] as long as it develops an adequate degree of reliability to operate in harsh space environments. Indeed, when used in space, both ARM cores and FPGA fabric embedded in SoCs are vulnerable to radiation-induced soft-errors [START_REF] Ebrahimi | Comprehensive Analysis of Sequential and Combinational Soft Errors in an Embedded Processor[END_REF][START_REF] Kastensmidt | Fault-Tolerance Techniques for SRAM-Based FPGAs[END_REF], which pose a greater reliability threat to SRAM-based FPGAs, such as those from Xilinx and Altera. In the latter FPGAs, the charged particles and outer radiation in general can alter the configuration information stored in SRAM-based memory cells, resulting in undesired logic functions implemented in the programmable logic and/or wrong inter-connections between the components. On the other hand Microsemi uses flash memory in its FPGA fabric, which is more resilient to radiation provoked soft-errors but allows for lower integration density, thus delivering more modest computation capabilities. Scrubbing is a classical method to protect the configuration memory in SRAM-based FPGAs. This technique consists in periodically checking the Error Correction Codes (ECCs) associated to the configuration information stored in the FPGA configuration memory and correct any errors that might have been occurred by rewriting the correct value, which is typically stored in an external rad-hard non-volatile flash memory. That said, Xilinx has released several generations of radiation-hardened FPGAs (e.g., Virtex-5QV [START_REF]Radiation-Hardened, Space Grade Virtex5-QV FPGA Data Sheet: DC and AC Switching Characteristics[END_REF]) and software tools for making designs fault-tolerant (e.g., Xilinx TMR Tool [15]) that are used in a number of space systems. The Xilinx roadmap includes the development of a radiation-tolerant SoC technology as well as the necessary software tools for creating fault-tolerant designs on it.

Current space instrument payload systems typically include either a flashbased FPGA (e.g., Microsemi ProASIC3 [START_REF] Inc | Radiation-Tolerant ProASIC3 Low Power Spaceflight Flash FP-GAs with Flash Freeze Technology[END_REF]) or a rad-hard SRAM-based FPGA (e.g., Xilinx Virtex5-QV) for implementing data acquisition, synchronization and processing, and an antifuse-based FPGA (e.g., Microsemi RTAX [START_REF] Inc | RTAX-S/SL and RTAX-DSP Radiation-Tolerant FPGAs[END_REF]) for implementing data communications with spacecraft main computer, internal bus handling, housekeeping data collection and management of the configuration of the SRAM-based FPGA. In applications that are critical for spacecraft mission, such as Guidance Navigation and Control (GNC), the antifuse FPGA is replaced by a rad-hard processor such as a BAE Systems RAD750 [START_REF]RAD750 Radiation-Hardened PowerPC Microprocessor[END_REF] or a Cobham Gaisler Leon3 [START_REF] Gaisler | LEON3-FT SPARC V8 Processor Data Sheet and User's Manual[END_REF]. A couple of recent NASA instruments that use this classic architecture are the ChemCAM on the Mars Curiosity rover [START_REF] Wiens | The ChemCam Instrument Suite on the Mars Science Laboratory (MSL) Rover: Body Unit and Combined System Tests[END_REF] and the Goddard Space Flight Center (GSFC) SpaceCube [START_REF] Petrick | Adapting the SpaceCube v2.0 Data Processing System for Mission-Unique Application Requirements[END_REF]. Hence, a SoC that includes these two components (processor + programmable logic) into a single chip is perfectly suited for space instrument payload systems.

Two are the reasons that have made us choose Xilinx Zynq SoC to prototype our APEX-SoC platform. First, Xilinx is one of the vendors with the most advanced SoC technology roadmap, which also addressed radiation-hardened FPGAs. Second and most important, NASA has recently approached Xilinx technology in the scope of its CubeSat Launch initiative (CSLI), as described in the paragraph below. Xilinx Zynq SoCs integrate a dual-core ARM Cortex-A9 centric Processing System (PS) and a 28 nm Xilinx 7-Series (Artix-7 or Kintex-7) Programmable Logic (PL) fabric. The chip includes abundant on-chip AXI ports with low power rails to communicate the PS with the PL, which results in substantially less power consumption, considerably higher bandwidth and lower latency.

The Xilinx Zynq SoC is in the heart of the Computer Space Processor (CSP) designed by the National Science Foundation (NSF) Center for Highperformance Reconfigurable Computing (CHREC) and licensed for fabrication to Space Micro Inc. [START_REF] Rudolph | CHREC Space Processor: a Multifaceted Hybrid Architecture for Space Computing[END_REF][20][START_REF] Mandl | Intelligent Payload Module Update[END_REF]. The CSP uses a combination of commercial and rad-hard components, where commercial devices perform critical computations and are supervised by the rad-hard devices (e.g., reset and watchdog circuits). This Zynq-based processor will be part of future NASA missions such as the Space test Program-Houston-ISS-5 SpaceCube experiment [START_REF] Wilson | CSP Hybrid Space Computing for STP-H5/ISEM on ISS[END_REF] and the Compact Radiation bElt Explorer (CeREs) heliophysics CubeSat [START_REF] Kanekal | CeREs: a Compact Radiation bElt Explorer[END_REF]. PlanetiQ Inc. will also integrate 3 CSPs on each of the 12 LEO weather satellites scheduled to be launched in 2017. In addition to these space missions, the CSP has been tested in neutron radiation and heavy-ion environment by Brigham Young University [START_REF] Wirthlin | Neutron Radiation Test Results of the Linux Operating System Executing within the CHREC Space Processor (CSP)[END_REF]. JPL, Xilinx and Swift LLC have also tested the Xilinx SoC part and other Xilinx 7-series FPGAs under heavy-ions radiation [START_REF] Amrbar | Heavy Ion Single Event Effect Measurements of Xilinx Zynq-7000 FPGA[END_REF][START_REF] Switft | Investigation of High Current Events in 28nm 7-series FPGAs[END_REF][START_REF] Koszek | Challenges in Assessing Single Event Upset Impact on Processor Systems[END_REF].

The APEX-SoC Platform and Infrastructure

The APEX-SoC platform is currently prototyped on a ZedBoard mini-ITX board, which is populated with a Xilinx Zynq 7Z100 SoC device. An FMC board containing an ADC is attached to the ZedBoard to deal with the analog electrical signals that are typically delivered by space science instruments. The APEX-SoC platform is also coupled with two external DDR memories to enable its use with instruments that generate large amounts of data: the PS-DDR is solely dedicated to the ARM processor in the Zynq, while the PL-DDR is used as scratchpad memory by the data processing modules implemented on the FPGA fabric and is also accessible by the ARM processor to retrieve the intermediate results computed by these. The last external component connected to the APEX-SoC is a SATA Solid State Device (SSD). This is used to temporarily store the (likely large amounts of) science results produced by the APEX-SoC until a downlink communication window with Earth is available, allowing for creating independent and stand-alone instruments avionics subsystems. The typical data-flow in the APEX-SoC is thus as follows: 1) the instrument data is acquired and processed by the FPGA logic, 2) the computed intermediate results by the FPGA logic are DMA-transferred to the DDR memory dedicated to the ARM processor for final processing, and 3) the final results are copied to the SSD prior to being downloaded to Earth.

The APEX-SoC provides support for integrating multiple identical data processing stages that can be used to process different science data in parallel to increase performance, or to detect computation errors by comparing their results when they process the same science data. This flexibility is needed when the requirements might change during the mission.

Fig. 1 shows a block diagram of the APEX-SoC architecture. The following subsections describe the major aspects related to this architecture as well as the main fault-tolerance mechanisms that are implemented on it.

ARM-centric Processing System

The ARM-centric PS includes all the peripherals that are typically required by flight science instruments, including: DMA support, GPIOs, Ethernet, SATA, interrupt controller and a memory-mapped register bank to exchange state and configuration data with the FPGA processing logic. As previously mentioned, process data are exchanged with the FPGA logic through the DMA-accessible PL-DDR memory. In order to speed-up the development of APEX-SoC-based instruments avionics, one of the ARM cores runs a standard Linux-based operating system, which provides Ethernet protocol to communicate with the spacecraft's main computer and a file system to ease the management of science results stored in the SSD. The second ARM core can be dedicated for software-based processing of instrument data. One scenario where software processing is convenient is when dealing with floating-point intensive algorithms, which can be easily computed using the NEON Floating Point Unit (FPU) [START_REF] Arm Ltd | Cortex-A9 NEON Media Processing Engine: Technical Reference Manual[END_REF] available in the ARM processor. A Real-Time Operating System (RTOS) can be deployed in this core Fig. 1: The APEX-SoC platform to use software multitasking to extend the hardware parallel processing carried out in the FPGA fabric while ensuring a sustainable use of CPU by all of the tasks [START_REF] Inc | Zynq All Programmable SoC Linux-FreeRTOS AMP Guide[END_REF].

Data-Flow Infrastructure

Each data processing module in the FPGA fabric is assigned a private data segment in the PL-DDR, with its size depending on the computing needs of that particular module. In order to exploit the full bandwidth delivered by the PL-DDR memory (6.4 GB/s) and to support the parallel/redundant execution of the hardware modules, the APEX-SoC implements eight 32-bit DDR access ports at 200 MHz using Xilinx-provided AXI-Stream Data Movers, which act as DMA controllers for the FPGA processing logic [START_REF] Inc | AXI DataMover User Guide[END_REF]. One of the DDR ports is dedicated to the instrument data acquisition logic (shown in blue color), another one is assigned to the ARM DMA, and the remaining six ports are connected to a crossbar that multiplexes them among the instrument data processing stages. The objective of this crossbar is thus to create as many communication channels as needed by the instrument-dependent modules using the physically available DDR ports. A data-flow controller drives the connections in the crossbar and schedules the PL-DDR accesses to maximize performance. For each data transfer, it specifies the memory address and size of the data segment to be read or written to the corresponding Data Mover. The data-flow controller is based on a tiny Xilinx 8-bit PicoBlaze processor [START_REF] Chapman | PicoBlaze for Spartan-6, Virtex-6, 7-Series, Zynq and UltraScale Devices[END_REF], which consumes only 26 LUTs in the Zynq FPGA fabric, and implements a collection of reusable assembler routines that provide the required flexibility to deal with a wide range of instruments. Most of the HDL code used to describe the APEX-SoC infrastructure is also parameterizable and can be easily customized to the needs of any instrument.

Fault-Tolerance Features

The temperature on the Zynq die is continuously monitored using an on-chip sensor (see XADC in Fig. 1) [START_REF]LogiCORE IP AXI XADC v1.00a Product Guide[END_REF] to identify and prevent overheat situations that could lead to the eventual destruction of the chip. Excessive noise situations in the power supply are also detected with this sensor. These may indicate that there is a problem with the voltage regulators, power lines in the PCB or even in the spacecraft power subsystem. Finally, the PL-DDR AXI Stream ports are continuously monitored to detect stuck-at situations and errors in memory data transfers. All storage resources in the APEX-SoC platform are protected with ECCs. The Xilinx ECC solution built in the silicon of the Zynq is used for the PS-DDR, whereas a custom ECC logic for the PL-DDR is implemented on the FPGA fabric. This ECC logic uses Hamming [START_REF]LogiCORE IP AXI XADC v1.00a Product Guide[END_REF][START_REF] Switft | Investigation of High Current Events in 28nm 7-series FPGAs[END_REF] codes to protect the data words transferred through each of the PL-DDR ports and is pipelined to maximize performance. It allows for detecting and automatically correcting single bit flips (e.g., radiation-induced SEUs) in a PL-DDR data word and detecting, but not correcting, double bit errors. Note that the possibility that multiple bit errors are accumulated in the same data word is small, as the ECC logic corrects every single bit flip that might have occurred in the short period of time data remains stored in the PL-DDR memory between consecutive write accesses. Data words affected by uncorrectable double bit errors can be either replaced by zeros or with the interpolated value of the two neighboring samples. All the finite state machines in the APEX-SoC are implemented using "one-hot" encoding, in such a way that radiation-provoked upsets in state flip-flops result in the state machine flow being redirected to an "illegal" state that signals the ARM processor the error situation. The correctness of the configuration data stored in the Zynq configuration memory is periodically checked by a Xilinx Single Event Mitigation (SEM) controller [START_REF]Soft Error Mitigation Controller[END_REF]. Single-bit upsets are automatically fixed by the Xilinx SEM, and in the event of a double bit upset, the ARM processor carries out a full reconfiguration of the FPGA fabric.

Reliability Mode

As previously introduced, the APEX-SoC permits to increase system reliability by using multiple identical data processing stages in an N-out-of-M scheme. The number of M redundant stages that can be implemented is only limited by the amount of FPGA resources available on the fabric and the energy budget, however Dual Modular Redundancy (DMR) or Triple Modular Redundancy (TMR) are typically used. In all cases, three redundant copies of the same science data are kept in the PL-DDR memory and replicated majority voters are connected both at the input and output of the M redundant processing stages as shown in Fig. 2. The input voters do not consider corrupted data that cannot be recovered using ECCs. When any of the output voters detect that all of its input results are different, a computation error is assumed and the processing of that science dataset is repeated. Computation errors can occur when radiation affects data registers and/or FPGA configuration [START_REF] Kastensmidt | Fault-Tolerance Techniques for SRAM-Based FPGAs[END_REF]. While upsets in the data registers cannot be detected by the Xilinx SEM controller, these are automatically corrected when reloading the data to process again. The Xilinx SEM is still needed to deal with the corrupted configuration bits, as described in section 3.3 The data-flow controller coordinates the access by the voters to the redundant data in the PL-DDR in a ping-pong fashion, so that the voted results do not overwrite the source data, in case the computation needs to be repeated.

Case-Study: APEX-SoC-based Controller of the JPL CIRIS Spectrometer

This section describes a proof-of-concept SoC implementation of a controller for the JPL CIRIS spectrometer using the APEX-SoC platform.

The JPL CIRIS Spectrometer

CIRIS is one of the new generation JPL instruments proposed to search for life indicators in icy moons, such as Europa [START_REF] Carlson | The Compositional InfraRed Interferometric Spectrometer (CIRIS) for Assessing the Habitability of Europa[END_REF]. It is based on the COTS instrument prototype described in [START_REF] Wadsworth | Rugged High-Speed Rotary Imaging Fourier Transform Spectrometer for Industrial Use[END_REF], and it a small, rugged and lightweight Fourier More importantly, as opposed to related instruments such as grating spectrometers (e.g., Galileo NIMS [START_REF] Carlson | Near-Infrared Mapping Spectrometer Experiment on Galileo[END_REF]), CIRIS has intrinsic immunity from radiationinduced noise, enabling it to perform mid-IR solar reflectance and thermal emission spectroscopy with limited interference from the radiation environment in Europa.

The major structural novelty introduced by CIRIS is the constant-velocity rotating refractor it uses to vary the optical path difference of the two rays in which incoming light is divided by a beam splitter at the entrance of the instrument (red and green rays in Fig. 3). The reflected rays in the rotating refractor recombine after travelling through the instrument, resulting in a fringe interference light pattern (interferogram) that is measured with a photo-detector (purple ray in Fig. 3). There are up to four regions over the course of a revolution of the refractor where the optical interference between the input light rays can be measured with a photo-detector. These regions are located at approximately 16 • arcs around the four positions where the refractor is parallel or perpendicular to the beam splitter. Note that the interferogram amplitude value is maximum in these four positions as all of the light rays travel the same distance along the spectrometer and recombine in phase at its output. This is why these positions are called Zero Path Difference (ZPD) positions. An optical incremental encoder mounted on the servomotor that drives the refractor's rotation is used on the ground prototype of CIRIS to identify these regions. As the CIRIS refractor performs 6.5 revolutions per second, each interferogram spans over a period of 13.6 ms every 24.8 ms. The optics and functioning of CIRIS result in an interferogram with the high-amplitude values assembled in a narrow central burst, and smallamplitude values spanning the vast majority of the tail positions and carrying the spectral resolution information (see fig. 4). The interferogram signal delivered by the photo-detector is conditioned, filtered and amplified to ±5V range prior to being digitized at 1 MSPS using the ADC available in the APEX-SoC. The interferogram samples are then processed via a Fast Fourier Transform (FFT) to produce a spectrum that illustrates the intensity of the wavelengths present in the light beam. This in turn permits to find out the chemical composition of the sample or body under study by looking at the absorption lines in the spectrum. However, spectral leakage (e.g., "picket-fence" effect) and noise are also present in the spectrum due to the limited discretization of the interferograms through time limited digital sampling, and need to be properly handled by the instrument electronics to produce meaningful results [START_REF] Saptari | Fourier Transform Spectroscopy Instrumentation Engineering[END_REF].

Fig. 3: CIRIS Spectrometer

Although radiation has small impact on the spectral content of CIRIS data, FTS data processing allows increasing the SNR of the instrument even further [START_REF] Yiu | Adaptive Controller for a Fourier Transform Spectrometer with Space Applications[END_REF]. As shown in fig. 4, the shape of the CIRIS interferogram allows the data processing for detecting (and removing) most of the radiation hits that induce large current pulses (i.e., significantly greater than the nominal value) in the instrument's photo-detectors. In fig. 4, note there are two radiation hits at -266 and -500 µs.

Fig. 4: CIRIS interferogram with radiation hits

At the moment there is a single photo-detector in the ground prototype of CIRIS, however the flight version of CIRIS will be equipped with an array of up to 25 photo-detectors to increase the instrument's spatial resolution and sensitivity in different IR bands. This will also increase the computation burden, as more interferograms will need to be processed within the same span of time (24.8 ms).

CIRIS Data Processing

The section describes the different processing stages that must be applied on the CIRIS interferogram data in order to produce meaningful spectroscopy results that can be interpreted by the scientists on Earth [START_REF] Saptari | Fourier Transform Spectroscopy Instrumentation Engineering[END_REF]. Fig. 5 shows a block diagram of these stages as well as their interfaces with the APEX-SoC infrastructure. In this figure, note the two superposed main blocks that represent the dual data processing solution adopted in the CIRIS APEX-SoC to increase the performance and reliability.

The first stage prepares the digitized interferogram samples for subsequent processing by selecting 8,192 samples centered around the ZPD positions. This is done to deal with any temporal shift that might have occurred while sampling the interferogram.

The second stage removes the DC offset in the ZPD aligned interferogram by subtracting its average value, which is computed using a Cumulative Moving Average (CMA).

The third stage implements a radiation hit filter to detect and remove the outlier in the interferogram provoked by radiation striking the CIRIS photodetector. The radiation pulses at the output of the CIRIS transconductance am-Fig. 5: CIRIS data processing block diagram plifier, with a bandwidth of 100 kHz, are about 10 µs full width at half maximum and easily recognized in the small-amplitude tail samples using statistics, namely the mean and variance (shown by triangles in Fig. 6) [START_REF] Yiu | Adaptive Controller for a Fourier Transform Spectrometer with Space Applications[END_REF]. Radiation hit samples are then replaced by zeros without modifying significantly the spectral content of the interferogram. This property comes from the fact that interferogram points outside the central burst mainly carry redundant resolution information, and hence, removing a few points out of 8,192 lead to indistinguishable changes in the spectrum. In effect, each of the interferogram samples contributes only in about 0.1% to the spectrum. Note here that the undetectable radiation hits that are at or below the un-irradiated noise level spread their energy over all wavelengths and therefore average to a constant DC offset in the spectrum, which is removed in the second stage. The mean and variance statistics are computed on the tail samples of the interferogram using the Knuth algorithm [START_REF] Knuth | The Art of Computer Programming[END_REF]. The fourth stage (STAT Inter.) computes the variance and performs a CMA on successive interferograms detected around the same ZPD positions with the objective of estimating and increasing the SNR by removing the effect of high frequency and random noise. As in the third stage, the Knuth algorithm is used to calculate these statistics.

The fifth stage apodizes the averaged interferograms at the edges of the sampled regions to minimize the effects of spectral leakage.

The sixth stage computes the FFT on the interferogram. In light of increasing spectral resolution, this stage adds 4,096 zeros to each of the tails of the interfer-ogram to obtain 8,192 additional interpolated spectrum points in-between the original nonzero-filled spectrum data, that is, 16,384 total spectrum points. This zero padding allows us to reduce the erroneous signal due to the "picket-fence" effect by up to 36% [START_REF] Gronholz | Understanding FT-IT Data Processing[END_REF].

The seventh stage relies on the Knuth algorithm to compute the variance and CMA on the spectrums resulting from the successive interferograms detected around the same ZPD positions.

The eighth and ninth stages are intended to correct the deviations provoked by CIRIS refractor's refractive index variations with wavelength and due to minor dissimilarities in the CIRIS optics between the four interferogram acquisition regions.

All of the processing stages described in this section are runtime configurable from Earth to adapt to potentially unexpected conditions when exploring distant planetary bodies. Some of the parameters that can be configured include: [START_REF] Blem | A Detailed Analysis of Contemporary ARM and x86 Architectures[END_REF] the method to detect the position of the ZPD sample (e.g., most-negative, mostpositive or most-magnitude) in the ZPD alignment core, (2) the number of trials to be averaged and the requirement to compute or not the variance in the STAT cores, (3) the apodizing function in the apodization core, and (4) the requirement for zero-filling (16,384 spectrum points) or not [START_REF] Iturbe | A Triple Core Lock-Step (TCLS) ARM Cortex-R5 Microprocessor for Safety-Critical and Ultra-Reliable Applications[END_REF]192 spectrum points) in the FFT core.

It is important to note here that the interferograms detected in the photodetector array that will be available in the flight version of CIRIS are independent of each other, and hence, the processing stages presented above are suitable for a parallel implementation on the APEX-SoC. Currently we simulate the photo-detector array by copying multiple times (to different PL-DDR memory segments) the same interferogram samples digitized by the single ADC in the system.

CIRIS Data Processing Integration into the APEX-SoC Infrastructure

All of the processing stages presented in section 4.2 are implemented on the FPGA fabric, except stages 7 and 8, which run as software routines in the ARM processor because they involve floating-point operations. Two instances of the whole IRIS data processing are integrated into the APEX-SoC infrastructure, as shown in Fig. 5. As previously mentioned, these can be used to boost performance or to improve reliability (i.e., DMR scheme), depending on the mission's requirement at each time. For this specific scenario, a crossbar with 13 communication channels is created and the associated assembler routines in the data-flow controller are appropriately tailored for the data transfers required by CIRIS processing stages.

Results

This section summarizes the implementation, performance and irradiation results we have collected so far in the APEX-SoC-based CIRIS controller.

Implementation

The amount and type of resources consumed by the APEX-SoC infrastructure and the CIRIS data acquisition and processing modules when implemented on a Zynq 7Z100 device are detailed in Table 1. The most important implementation aspect to note here is the spatial isolation of the CIRIS modules within the FPGA fabric, which has been carried out following the Xilinx Isolation Design Flow (IDF) [START_REF] Hallett | Isolation Design Flow for Xilinx 7 Series FPGAs or Zynq-7000 AP SoCs (ISE Tools)[END_REF]. This permits to increase the availability of the system by preventing the situation where a single charged particle corrupts multiple processing stages in the FPGA fabric. As shown in Fig. 7, the crossbar and the data-flow controller are mapped in-between the two processing stages forming a fence. It is also important to note the small footprint of the voters compared to the processing stages, in the range of hundreds of LUTs and flip-flops, which minimizes the chances of being corrupted by radiation. The power consumption reported by Xilinx Vivado design tool for the whole APEX-SoC-based CIRIS controller is approximately 5 W, which is about 3 W less than that of an equivalent board based controller.

The APEX-SoC uses up to 256 MB in the PL-DDR memory (approximately 25% of the total DDR memory capacity) to process 25 interferograms simultaneously. The memory is arranged into several data segments across different categories, each containing a given type of data (e.g., raw interferogram, interferogram mean/variance, spectrum amplitude mean/variance or spectrum phase mean/variance) related to the information detected by a given photo-detector when the rotating refractor was in a given position. Besides, as explained in section 3.4, each data is stored three times in the PL-DDR memory in different TMR data pools to increase reliability, and each TMR pool is itself replicated two times (A and B) to allow for data re-processing, if needed.

Performance

Table 2 shows the performance results measured when the FPGA processing logic is clocked at 200 MHz and both DDR memories and ARM Cortex-A9 processor run at 800 MHz. As shown in Fig. 8, the parallelism provided by the dual data processing channels in the APEX-SoC and the efficient PL-DDR memory access schedule allow for processing two interferograms every 460 µs, with a latency allocated to processing and the long latency introduced by the radiation hit filter and the FFT computation is hidden by overlapping parallel processing and PL-DDR data transfers. As a result, the APEX-SoC almost quadruples the processing requirements of the flight CIRIS spectrometer as it is able to process about a hundred interferograms within the time span the refractor in the instrument takes to get between consecutive ZPD positions (24.8 ms). On the other hand, when using the two redundant data processing channels to increase reliability (i.e., DMR scheme), the APEX-SoC-based CIRIS controller is still able to fulfill the processing requirement for the next-generation of CIRIS, requiring up to 900 µs to process each interferogram.

CIRIS Processing

Receive

Robustness Against Radiation

A radiation test was conducted at JPL using a 60Co γ-ray source (1 rad/sec) directed toward the CIRIS photo-detector operating at 77K to reduce detector noise below the radiation hit pulses. The hit rate in this test was approximately 3,400 hits per second, exceeding what is expected in the Europa mission by at least a factor of three. Fig. 9 shows the obtained results, where the blue line represents the measured values without radiation, the purple line the values measured with radiation and the yellow line represents the values measured when the radiation hit filter was enabled. As shown in Fig. 9a, the radiation hit filter stage reduces the distortion of the line shape and spectrum provoked by radiation while keeping all other spectral components unaltered. Note in this figure that the high artificial "emission" peaks on the spectrum are coming from electrical noise generated by the vacuum pumps in the laboratory. These results are of utmost importance towards building an instrument that could cope with Europa-like radiation, which indeed deteriorated the spectroscopy data collected by NASA's previous generation NIMS spectrometer aboard the Galileo spacecraft more than a decade ago [START_REF] Carlson | Near-Infrared Mapping Spectrometer Experiment on Galileo[END_REF]. In addition, as shown in Fig. 9b, the adopted mitigation solution increases the instrument SNR by eliminating the noise due to radiation hit pulses.

We have not conducted any specific experiment to test the implemented faulttolerance features yet, as these are well known and proven to be effective. Plans are to port the APEX-SoC-based CIRIS instrument described in this chapter to a radiation-hardened Xilinx SoC as soon as this technology is available and test the design in a simulated Europa-like thermal and radiation environment.

Conclusions and Future Work

This chapter has presented an ongoing research conducted by NASA's Jet Propulsion Laboratory (JPL) and ARM to develop a SoC platform (APEX-SoC) to power instruments avionics in future space exploration missions. This platform reduces significantly the size and power consumption of the instrument avionics as most of the electronics required for science processing of the instrument data are fitted in a single chip. At the moment this platform is prototyped using a commercial Xilinx Zynq SoC, where a number of fault-tolerance mechanisms have been implemented. The expectation is to port this design to a radiationtolerant SoC part that might be available in the near future. The chapter has presented a case-study where the APEX-SoC prototype is used to process data delivered by a JPL spectrometer (CIRIS). Finally, the chapter has discussed the implications of using SoC technology in future space missions.

Future work at ARM will focus on making the processor more resilient to radiation. Namely, a thorough study of the ARM Cortex-R5 microarchitecture will be conducted to identify the parts that are more vulnerable and to choose the most suitable fault-tolerance techniques to be used in each of these parts without compromising the area and power consumption efficiency of the ARM At the processor architecture level, a fail-operational Triple Cortex-R5 Core Lock-Step (TCLS) processor will be developed [START_REF] Iturbe | A Triple Core Lock-Step (TCLS) ARM Cortex-R5 Microprocessor for Safety-Critical and Ultra-Reliable Applications[END_REF]. JPL will look forward using the APEX-SoC platform with other space instruments.

Besides the research described in this chapter to design the next-generation space instruments avionics, JPL is also working in collaboration with the Goddard Space Flight Center (GSFC) and the Air Force Research Laboratory (AFRL) on designing a next-generation high-performance spaceflight processor based on a dual quad-core ARM Cortex-A53 [START_REF]Statement of Work (SOW) For the Development of the High Performance Space Computing (HPSC) Processor[END_REF].

Fig. 2 :

 2 Fig. 2: DMR scheme implemented in the APEX-SoC-based CIRIS controller

  Interferograms at Position 1 at Line 1 DŝƚŝŐĂƚŝŽŶ ŽĨ ϭ ƌĂĚͬƐĞĐ ZĂĚŝĂƚŝŽŶ

Fig. 6 :

 6 Fig. 6: CIRIS interferogram radiation effects and mitigation

FenceFig. 7 :

 7 Fig. 7: APEX-SoC-based CIRIS controller floor-planning

Fig. 8 :

 8 Fig. 8: Timing diagram and PL-DDR schedule (STATi=STAT Interferogram, STATa=STAT Amplitude STATp=STAT Phase)

  Radiation Spectra at position 1 line 1 across 128 trials beginning at line 1 with windows of 33 points (b) SNR improvement

Fig. 9 :

 9 Fig. 9: Radiation mitigation in APEX-SoC-based CIRIS data processing

Table 1 :

 1 Resources consumed in a Xilinx Zynq 7Z100 SoC

	Component	LUTs	Flip-flops	DSP48s BRAM36s
	Data Acquisition	347	267	N/A	N/A
	RHBD Features	6,194	3,818	N/A	14.5
	Data Processing	61,874	46,152	414	155
	Infrastructure	61,773	49,631	N/A	395
	Total	130,148 (47%) 99,868 (18%) 414 (20%) 564.5 (75%)

Data Process Data Deliver Data

	ZPD Alignment	100 µs	N/A	41 µs
	DC Removal	N/A	281 µs	N/A
	Radiation Hit Filter	N/A	320 µs	N/A
	STAT Interferogram	229 µs	∼0 µs	188 µs / 41 µs
	Apodization	47 µs	∼0 µs	41 µs
	FFT	41 µs	165 µs	43 µs
	STAT Spectrum	229 µs	∼0 µs	188 µs

Table 2 :

 2 Performance results

Acknowledgment

The research described in this chapter was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration (NASA). Xabier Iturbe is funded by the European Commission's FP7 Marie-Curie International Outgoing Fellowship Program with "Project No. 627579".


Mark Strugatsky

Kira Seleznyova

Sergey Yagupov

Alexey Drovosekov

Janis Kliava

Nature of magnetocrystalline anisotropy in the basal plane of iron borate

Keywords:

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

Introduction

From the standpoint of magnetic properties, iron borate FeBO 3 is an easy plane antiferromagnet, possessing a weak in plane fer romagnetic moment. Iron borate has rhombohedral calcite struc ture of point group symmetry D 3d and the space group D 6 3d [START_REF] Diehl | Current Topics in Materials Science[END_REF]. The effective basal (hexagonal) magnetocrystalline anisotropy energy for FeBO 3 can be expressed as follows [START_REF] Strugatsky | [END_REF]:

E A e eff cos 6u ð1Þ
where u is the azimuthal angle of the antiferromagnetic vector and e eff is the effective basal anisotropy constant,

e eff e FeBO 3 þ 1 4

d 2 FeBO 3 a eff ; ð2Þ
d FeBO 3 and e FeBO 3 being the basal magnetocrystalline anisotropy con stants and a eff being the effective uniaxial magnetocrystalline aniso tropy constant,

a eff a FeBO 3 þ D 2 FeBO 3 E FeBO 3 : ð3Þ
In the latter equation a FeBO 3 ; D FeBO 3 and E FeBO 3 are, respectively, the uniaxial magnetocrystalline anisotropy, Dzyaloshinskii Moriya and exchange constants for FeBO 3 .

As far as for Fe 3+ (3d 5 electronic configuration) the orbital moment equals zero, the exchange energy in a good approximation is isotropic [START_REF] Kanamori | Magnetism, A Treatise on Modern Theory and Materials[END_REF], so that the magnetocrystalline anisotropy energy of FeBO 3 includes only crystal field (cf) and dipole dipole (dip) terms:

a FeBO 3 a cf þ a dip d FeBO 3 d cf þ d dip e FeBO 3 e cf þ e dip ð4Þ
The crystal field contributions to the abovementioned constants have been calculated by Seleznev in the mean field approximation [START_REF] Seleznev | [END_REF]. Indeed, the effective exchange field H E in iron borate is in order of 10 3 kOe, thus, Fe 3+ ions experience strong exchange coupling [START_REF] Baryahtar | Abstracts of 19th Soviet Symp. on Low Temperature Physics[END_REF]. The crystal field contribution to the anisotropic part of the magne tocrystalline anisotropy energy can be calculated in perturbation theory using the spin Hamiltonian for isolated Fe 3+ ion in a dia magnetic crystal isomorphous to iron borate, e.g., gallium borate [START_REF] Yagupov | [END_REF]7]. Thus, a cf ; d cf and e cf can be expressed through the parame ters of the spin Hamiltonian. The corresponding expressions will be given and the parameters will be specified in the next Section.

The dipole dipole contributions to these constants are usually calculated using the lattice sum method. The value of a dip 3:82 10 5 Jm 3 at 0 K for FeBO 3 has been obtained previously [START_REF] Rudenko | Abstracts of Soviet Symp. on Physics of Magnetic Phenomena[END_REF]. a FeBO 3 has been determined by antiferromagnetic resonance ⇑ Corresponding author. E-mail address: janis.kliava@u-bordeaux.fr (J. Kliava).

1

(AFMR) in a wide temperature range by Velikov et al. [START_REF] Velikov | [END_REF]; its value extrapolated to 0 K is a exp Earlier, the occurrence of the dipole dipole contribution to hexagonal basal anisotropy in iron borate had been ruled out on the grounds of symmetry. Indeed, the dipole dipole interaction energy is usually considered for ''point dipoles" having a negligible size, in which case only the uniaxial anisotropy is accounted for. Meanwhile, more sophisticated anisotropies, in particular, the hexagonal magnetocrystalline anisotropy can be described by higher than first order terms in the expansion of the dipole dipole interaction energy in a Taylor series in the small parameter dipole size/interdipole distance. Thus, taking into consideration ''extended dipoles", having non negligible size, opens the possibil ity to reasonably account for the dipole dipole contribution to the basal anisotropy constants and, subsequently, to estimate effective dipole dimensions in iron borate.

With this aim in mind, we have developed a theoretical descrip tion of two models of the extended magnetic dipole: (1) an assem bly of two fictitious ''magnetic charges" AEq a distance d apart and (2) an Ampérian current, i.e., a circular current loop of a radius R delimiting an area S pR 2 . Previously, we have obtained and dis cussed the expressions of vector potentials and magnetic fields produced by these two models [10]. The purpose of the present work has been to calculate the dipole dipole contribution to the basal magnetocrystalline anisotropy constants of iron borate and to evaluate the size of dipoles in FeBO 3 in the framework of the suggested models.

Crystal field contribution to magnetocrystalline anisotropy

In order to calculate the crystal field contribution we have to consider a Fe 3+ ion in a diamagnetic crystal isomorphous with FeBO 3 . The conventional spin Hamiltonian in this case is [11,12]:

H gbH Á S þ 1 3 DO 0 2 1 180 ða FÞO 0 4 þ 2 p 9 aðÇO 3 4 sin 3a þ O 3 4 cos 3aÞ ð
5Þ
where g is close to the free electron g factor g e 2:0023, b is the Bohr magneton, H is the magnetizing field, S 5 2 is the electron spin of Fe 3+ , D is the second order axial fine structure constant, a and F are, respectively, the fourth order cubic and axial fine structure constants and O 0 2 ; O 0 4 ; O 3 4 and O 3 4 are extended Stevens operators, as defined in the textbook by Al'tshuler and Kozyrev [START_REF] Al'tshuler | EPR in Compounds of Transition elements[END_REF]. The Ç signs refer to two non equivalent iron sites with local magnetic axes rotated through the angle Ça about the C 3 axis, see Fig. 1.

The spin Hamiltonian parameters of Fe 3+ in diamagnetic GaBO 3 have been determined previously by Lukin et al. [11] and recently specified by Seleznyova et al. [12].

Seleznev has obtained the expressions for a cf ; d cf and e cf using the mean field approximation [START_REF] Seleznev | [END_REF]:

a cf 2N½D tðxÞ 1 2 ða FÞrðxÞ d cf 2 p 3 N a rðxÞ cos 3a e cf 0 ð6Þ
where N is the number of Fe 3+ ions per unit volume, N 2:236 10 28 m 3 for FeBO 3 and tðxÞ and rðxÞ are the following functions:

tðxÞ ð8ch 1 2 xþ2ch 3 2 x 10ch 5 2 xÞsh 1 2
x sh3x and rðxÞ 5

ð2ch 1 2 x 3ch 3 2 xþch 5 2 xÞsh 1 2 x sh3x : ð7Þ
where x gbH E =2kT, H E is an effective exchange field, vide infra, k is the Boltzmann constant and T is the absolute temperature.

As one can see from Eq. ( 6), the crystal field gives no contribu tion to e FeBO 3 .

With the spin Hamiltonian constants D 0:1032; a 0:0158; F 0:0368 cm 1 and a 36 o specified recently [12],

we obtain the following crystal field contributions to the magne tocrystalline anisotropy constants for FeBO 3 at 0 K:

a cf 4:82 Á 10 5 Jm 3 d cf 2:55 Á 10 3 Jm 3 e cf 0 : ð8Þ
3. Dipole-dipole contribution to magnetocrystalline anisotropy

Model of a pair of fictitious magnetic charges

Here we consider the interaction energy E between two paral lel/antiparallel identical dipoles for different models. Fig. 2 shows a system of two interacting dipoles implemented as a pair of Fig. 1. Two non-equivalent sites of Fe 3+ . ðx; y; zÞ is a Cartesian coordinate system with xkC2; ykm; zkC3 where C2; m and C3 are, respectively, a two-fold axis, a symmetry plane and the three-fold axis [START_REF] Diehl | Current Topics in Materials Science[END_REF]. The z-axis is perpendicular to the plane of the figure and points towards the reader. The full and empty circles represent ions located above and below this plane, respectively. ''magnetic charges" AEq spaced a distance d apart. For the dipole dipole energy one gets:

E AE l 0 4p m 2 d 2 2 r 1 r mp 1 r pm ð9Þ
where l 0 is the permeability of vacuum, m qd qde m is the mag netic moment, defined by analogy with electrostatics and directed along the unit vector e m ðsin#cosu; sin#sinu; cos#Þ;

r ðr 2 x þ r 2 y þ r 2 z Þ 1=2 is the distance between the centers of the dipoles, r mp r 2 þ d 2 þ 2rde r e m q and r pm r 2 þ d 2 2rde r e m
q where e r is the unit vector in the direction of r, see Fig. 2 for the notation. The choice of the AE signs refers to parallel and antiparallel dipoles, respectively.

Introducing e d=r in Eq. ( 9) one gets:

E AE l 0 4p  m 2 d 2 r 2 1 1 þ e 2 þ 2ee r Á e m p 1 1 þ e 2 2ee r Á e m p ! ð10Þ
For our purpose we need only an approximate expression of E for e ( 1. An expansion in the Taylor series up to the fourth order yields:

E Ç2 l 0 4p m 2 r 3 ðP 2 þ P 4 e 2 þ P 6 e 4 Þ ð 11Þ
where P n are Legendre polynomials [START_REF] Korn | Mathematical Handbook for Scientists and Engineers[END_REF] of the scalar product e r e m , and the and + signs correspond to parallel and antiparallel dipoles, respectively.

Model of a circular current loop

Next, we calculate the interaction energy between two identical and parallel circular current loops (Ampérian currents) of the same radii R and areae S pR 2 , carrying a current I, see Fig. 3. By defini tion, the magnetic moment of a loop is m SIe m pR 2 Ie m .

For definiteness, we choose the loops centered at the space ori gin O p and at an arbitrary point O s as the primary and secondary loops, respectively.

The dipole dipole interaction energy in this model is straight related to the mutual inductance M of the loops:

E MI 2 : ð12Þ
Here the sign in the right hand member occurs because the interaction energy between two coaxial (attracting) currents is negative while M in this case is positive. By definition,

M U I ð13Þ
where U, the magnetic flux induced by the current in the primary loop and passing through the secondary loop, can be calculated as follows:

U

I ls I lp dA Á dl s : ð14Þ
Here l p and l s are perimeters of the primary and secondary loops, dA is a differential element of the vector potential at a point M s on the secondary loop, produced by the primary loop: dA l 0 4p

I M p M s dl p ; ð15Þ
M p M s is the distance between two points on the primary and sec ondary loops, and dl p and dl s are differential elements of the corre sponding loops.

In order to evaluate the closed curve integrals in Eq. ( 14) we express M p M s , dl p and dl s as follows:

jM p M s j 2 r 2 þ 2R 2 ½1 sinð/ sÞ
þ2Rr x ½sin uðcos / þ sin sÞ þ cos u cos #ðsin / cos sÞ þ2Rr y ½sin u cos #ðsin / cos sÞ cos uðcos / þ sin sÞ 2Rr z sin #ðsin / cos sÞ ; In Eqs. ( 16) to (18) s and / are polar angles of arbitrary points of the primary and secondary loops, respectively. Putting these expressions in Eq. ( 14), for the interaction energy, cf. Eq. ( 12), we get:

dl p Rðcos u cos # sin s þ sin u cos sÞds Rð sin u cos # sin s þ cos u cos sÞds R sin # sin s ds 0 B @ 1 C A ð16Þ
and dl s Rðsin u sin / cos # cos u cos /Þd/ Rðsin u cos # cos / þ cos u sin /Þd/ R sin # cos / d/ 0 B @ 1 C A: ð17Þ
E l 0 4p I 2 Z 2p 0 Z 2p 0 dl p Á dl s M p M s ð19Þ
where the integrations are over s and /.

Since, as in the previous case, we need only an approximate expression of E, we can first expand the integrand in Eq. ( 19) in a Taylor series in the small parameter e R=r and then integrate the result. In terms of the Legendre polynomials, up to the fourth order we get:

E Ç2 l 0 4p m 2 r 3 P 2 3P 4 e 2 þ 75 8 P 6 e 4 ! ð20Þ
where the and + signs correspond to parallel and antiparallel dipoles, respectively.

Calculation of the dipole-dipole interaction energy in FeBO 3

In order to calculate the dipole dipole contribution to the mag netocrystalline anisotropy constants for FeBO 3 , we have put for ward a computer code implementing the lattice sum method. We have chosen to do the summation in the volume of a rhombo hedron congruent to the primitive rhombohedron shown in Fig. 4.

The axes of the rhombohedral coordinate system x 0 ; y 0 ; z 0 coin cide with the edges of the rhombohedron, see Fig. 4. In transform ing the radius vector from the Cartesian to the rhombohedral system, we express the coordinates of iron sites through the edge length l of the rhombohedron: x 0 m l, y 0 n l and z 0 k l, where m; n; k are integers numbering the sites along the corresponding axes. The radius vector in the new coordinate system is r l

ðz 0 x 0 Þ cos c k 0 i ðx 0 2y 0 þ z 0 Þ cos b i 0 j ðx 0 þ y 0 þ z 0 Þ cos c k 0 k 0 B B B @ 1 C C C A ð21Þ
where cos c k 0 i 1 cos b 2 q , cos b i 0 j 1 cos b 6 q and cos c k 0 k 1þ2 cos b 3
q are the cosines of angles between the corresponding axes of two coor dinate systems. Substituting Eq. ( 21) in the expressions for the dipole dipole interaction energy, Eqs. ( 11) and (20) for the models of a pair of magnetic charges and a circular current loop, respec tively, we express this energy in terms of integers numbering the iron sites along the edges of the rhombohedron.

The calculation of the dipole dipole energy density now is reduced to computing the following sum:

E dip 1 2 N X m;n;k ð 1Þ mþnþk Eðm; n; kÞ ð 22Þ
where the factor ð 1Þ mþnþk takes into account antiferromagnetic ordering and Eðm; n; kÞ is the dipole dipole interaction energy between ions at the origin (numbered 0; 0; 0) and at a site num bered m; n; k. Henceforth, the magnetic dipole moment at T 0 K will be expressed as

m gbS ð23Þ
where g; b and S have the same meanings as in Eq. ( 5).

The dipole dipole contributions at 0 K, together with those of the crystal field, described in Section II, are listed in Table 1. One can see that the models of a pair of magnetic charges and a circular current loop result in substantially different expressions for the dipole dipole energies.

In order to get the dipole dipole contributions to the magne tocrystalline anisotropy constants at different temperatures, these contribution at 0 K should be multiplied by ðM T =M 0 Þ 2 , where M T is the sublattice magnetization at the temperature T. ðM T =M 0 Þ for FeBO 3 have been tabulated [START_REF] Seleznev | [END_REF]. The temperature dependences of the crystal field contributions are given in Eq. ( 6).

Comparison with experiment

From the AFMR experiments the effective basal magnetocrys talline anisotropy constant e exp eff , cf. Eq. ( 2), can be determined [START_REF] Doroshev | [END_REF]. The EMR studies have been carried out at 77 K with a laboratory developed spectrometer at microwave frequencies m from 15 to 36 GHz and magnetizing field H up to 10 kOe applied in the basal plane of the crystal. FeBO 3 crystal has been synthesized Table 1 Crystal field and dipole-dipole contributions to the magnetocrystalline anisotropy constants of iron borate at 0 K. by solution in the melt technique and had the shape of a thin hexagonal plate [START_REF] Yagupov | [END_REF].

For H applied in the basal plane of the crystal, the low frequency AFMR mode for FeBO 3 is described by the following expression [START_REF] Doroshev | [END_REF]16]:

m c½HðH þ H D Þ þ 36H E H hex cos 6u þ H 2 D 1=2 ð24Þ
where c is the gyromagnetic ratio for g 2:0; H D and H E , H E H hex and This value of H D confirms the previous finding [START_REF] Velikov | [END_REF], see Table 2. Using this value, we have determined the product H E H hex , see Eq.

(24), from AFMR measurements at a fixed m, carried out by rotating H in the basal plane of the crystal.

Fig. 6 shows the angular dependence of the quantity

ðm=cÞ 2 HðH þ H D Þ.
One can see that this dependence can be accounted for by a superposition of hexagonal and uniaxial magnetocrystalline aniso tropies in the basal plane. The occurrence of uniaxial magnetocrys talline anisotropy in this case can be due to a slight deviation of H from the basal plane or caused by mechanical stresses [17]. There fore, the angular dependence of the resonance field has been fitted to by the following expression, cf. Eq. ( 24):

m c 2 HðH þ H D Þ H 2 D þ 36H E H hex cos 6u þ p cos 2u ð25Þ
where H 2 D ; H E H hex and p are fitting parameters, the term in p describing the uniaxial component. The best fit parameters

H 2
D and H E H hex are listed in Table 2 and p 1:13 AE 0:05kOe 2 .

From H E H hex and H E we get H hex 0:913 10 5 kOe; then, using the definition of H hex given in Table 2, we get e exp eff 0:936 Jm 3 . The H hex value corroborates that earlier reported by Doroshev et al. [START_REF] Doroshev | [END_REF]; note that the negative sign of the latter value is due to a dif ferent definition of the angle u. Finally, from the results described above we estimate the dipole size. In the following we assume that e eff e exp eff , see Eq. ( 2) for e eff . Substituting e exp eff in Eq. ( 2) and taking into account Eqs. ( 3) and ( 4), we get:

e exp eff e cf þ e dip þ 1 4
d cf þ d dip À Á 2 a FeBO 3 þ D 2 FeBO 3 E FeBO 3 : ð26Þ
In what follows, for a FeBO 3 we use the experimental value a exp FeBO 3 3:2 10 5 Jm 3 determined by AFMR at 77 K [START_REF] Seleznev | [END_REF][START_REF] Velikov | [END_REF].

D FeBO 3 and E FeBO 3 at 77 K have been calculated from the experimen tal values of corresponding effective fields and sublattice magneti zation, see Table 2. Substituting in Eq. ( 26) the abovementioned values as well as those given in Table 1, for the model of a pair of fictitious magnetic charges we get:

853 d 4 þ 135 d 2 þ 3:2 0:936 ð27Þ
Obviously, this equation can have only complex solutions; therefore, this model is not applicable in the actual case. In con trast, for the model of a circular current loop, with an analogous substitution, we get: 7:47 Á 10 3 R 4 405:4R 2 þ 3:2 0:936; ð28Þ yielding two positive solutions: R 1 0:2189 and R 2 0:0797Å: In order to assess the plausibility of such values, they should be compared with the ionic radius R i of Fe 3+ ; indeed, we can reason ably infer that the effective size of a dipole should be of the same order of magnitude as the size of the physical object producing the corresponding dipole moment. For high spin Fe 3+ in sixfold oxygen coordination R i 0:645 Å [18]; therefore the R 1 value seems to be a more realistic estimate than R 2 , the latter value being an order of magnitude smaller than R i .

The dipole dipole contributions to the magnetocrystalline ani sotropy constants of FeBO 3 at 77 K calculated with the R 1 value are: e dip 12:2 Jm 3 and d dip 7:57 Á 10 3 Jm 3 : ð29Þ

Conclusions

Possible contributions to the basal magnetocrystalline aniso tropy of iron borate, namely, crystal field and dipole dipole contri butions have been considered in detail. The former contribution has been calculated using the spin Hamiltonian parameters for iso lated Fe 3+ ions in (diamagnetic) gallium borate. The latter contribu tion has been evaluated assuming that the ratio dipole size/ interdipole distance is non negligible, i.e., that we are dealing with extended dipoles. The dipole dipole interaction energy has been calculated for two extended dipole models, viz., a pair of magnetic charges and a circular current loop. The dipole dipole contribution has been calculated by the lattice sum method.

In order to determine the basal magnetocrystalline anisotropy constants of iron borate, we have carried out AFMR studies at 77 K. A comparison between the experimental and calculated val ues of this constant has shown that the model of a pair of magnetic charges fails at explaining the experimental results. In contrast, the model of circular current loop provides consistent evidence in sup port of the dipole dipole contribution to the basal magnetocrys talline anisotropy of iron borate and, incidentally, yields a more or less realistic estimate of the size of the magnetic dipole associ ated with Fe 3+ ion.

In spite of the fact, that the modeling considering extended dipoles, put forward in this work provides new insight in the nat ure of the basal magnetocrystalline anisotropy of iron borate, it is certainly oversimplified. More sophisticated (ab initio) calculations are necessary in order to get further insight in the spatial distribu tion of the magnetic field produced by paramagnetic ions at short and intermediate distances.

Fig. 3 .

 3 Fig. 3. System of two interacting Ampérian currents.

Fig. 2 .

 2 Fig. 2. System of two interacting dipoles in the model of a pair of magnetic charges.

Fig. 4 .

 4 Fig. 4. A primitive rhombohedron with the edge length l 3:6 Å and the apex angle b 79:9 used to calculate the dipole-dipole energy.

4 À9:19 Á 10 5 R 4 =l 4 2 H 2 D 4 .À0: 9 Á

 442249 Table Parameters in Eq. (24) and their experimental values for FeBO 3 . 25 ± 0.25 kOe 2 [9] 6.543 ± 0.032 kOe 2 [this work] Anisotropic energy gap H E H hex 1 10 5 kOe [15] 0:913 Á 10 5 kOe [this work] a For FeBO 3 M 0 520G and M 77 % 512:9G [4]. b Value at 0 K.

H 2 D 2 D 6 :

 226 are, respectively, effective Dzyaloshinskii Moriya and exchange fields, anisotropic and isotropic energy gaps, H hex being the effective field of basal magnetocrystalline anisotropy, and u is the angle between H and the x axis. As far as different authors use different definitions of H D , H E and H hex , the definitions used in this work are given in Table 2 together with their experimental values. The field sweep EMR spectra recorded at different m allow obtaining the relation between m and H. The experimental results shown in Fig. 5 can be convincingly fitted to Eq. (24), resulting in the following best fit parameters: H D 98:68 AE 0:52kOe and 36H E H hex cos 6u þ H 46 AE 0:34kOe 2 :

Fig. 6 .

 6 Fig. 6. Dependence of ð m c Þ 2 À HðH þ HDÞ on the angle u. The curve is a fitting according to Eq. (25).

Fig. 5 .

 5 Fig. 5. Dependence of the AFMR frequency on the magnetizing field for FeBO 3 crystal. The dashed curve is the best fit according to Eq. (24).

Acknowledgments

This work was partially supported by the Russian Foundation for Basic Research and the Ministry of Education, Science and Youth of the Republic of Crimea in the framework of scientific pro ject Grant no. 16 42 910593 ''p_a" and by the V.I. Vernadsky Crim ean Federal University Development Program for 2015 2024.


Modeling and algorithm on a class of mechanical systems with unilateral constraints

Keywords: cable -structure, artificial NNW, linear complementarity problem

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

Introduction

The analysis of cable systems in the geometrically nonlinear range had been the object of a number of studies, for example, Greenberg [START_REF] Greenberg | Inelastic analysis of suspension roof-structures[END_REF], Jonatowski and Birnstiel [START_REF] Jonatowski | Inelastic stiffened suspension space structures[END_REF], Murray and Willems [START_REF] Murray | Analysis of ineladtic suspension structures[END_REF] and Tene and Epstein [START_REF] Tene | Nonlinear space trusses with tension-only members[END_REF]. But the behavior of individual cables becomes highly nonlinear if slackening and plastic tensile strains intervene. This behavior of cable can be called physical nonlinearity. Because of this physical nonlinearity the cable-structure belongs to a class of mechanical systems with unilateral constraint.

To avoid loosening phenomenon of the cables for practical constructions, pre-tension of cables is always used in design. However, for some cable-structures loosening phenomenon occurs still under dynamical load. As a typical example, in 1963, the roof of station hall Villars was damaged under wind load [START_REF] Jawerth | Ein Beitrag zur Frage der Eigenschwingungen, windanfachenden Kraefte und aerodynamischen Stabilitaet bei haengenden Daechern[END_REF]. The model experiment has shown that when an unsymmetrical load acts on this structure, some cable members may be loosening and the roof may be floating. Therefore description with consideration of only geometric nonlinearity for cable-structure can be not realistic. So by modeling structures with cables, the behavior of unilateral constraint must be considered.

The consideration of nonsmooth behavior for the cable-structures may date back to the 1970s. In 1975, Maier et al. [START_REF] Maier | Energy approach inelastic cable-structure analysis[END_REF] had dealt with energy approach to inelastic cable-structure analysis, in which they considered geometric nonlinearity and physical nonlinearity of cable members. They had adopted piecewise linear constitutive idealizations, particularly in limit and elastoplastic-work-hardening analysis. Their aim was to alleviate difficulties by using of smooth models [START_REF] Maier | Energy approach inelastic cable-structure analysis[END_REF]. In 1976, Panagitopoulos had dealt with the inelastic, stressunilateral analysis of cable-structures undergoing large displacements [START_REF] Panagiotopoulos | A variational inequality approach to the inelastic stress-unilateral analysis of cable-structures[END_REF]. The response of the cable-structures to load-and initial strain-increments is described by a set of equations and inequalities, which are formulated compactly as variational inequalities. Two dual extremum principles with inequalities as subsidiary conditions are derived. The numerical calculation is performed iteratively by using the decomposition techniques of the multilevel optimization and applying to the resulting substructures the algorithms of non-linear optimization. Based on the above mentioned works, analysis of cable-structures may be included in the complementarity problems. But at that time that study was on static or quasi-static problems.

The mechanics with unilateral constraint as an important subclass of non smooth mechanics develops rapidly since 35 years [START_REF] Brogliato | Nonsmooth mechanics[END_REF]. This field consists of two parts in the area of variational inequalities in connection with convex energy functions [START_REF] Monteiro Marques | Differential inclusions in nonsmooth mechanical problems[END_REF][START_REF] Moreau | Application of convex analysis to some problems of dry friction[END_REF][START_REF] Panagiotopoulos | Inequality problems in mechanics and applications[END_REF][START_REF] Rockafellar | Convex analysis[END_REF] and the area of hemivariational inequalities in connection with non-convex energy function, which is only 15 years old [START_REF] Panagiotopoulos | Hemivariational inequalities[END_REF]. Most nonsmooth mechanical problems possess non-convex features. In [START_REF] Panagiotopoulos | Inequality problems in mechanics and applications[END_REF] Panagiotopoulos considered inequality problems, which led him consistently to the development of hemivariational inequalities [START_REF] Panagiotopoulos | Hemivariational inequalities[END_REF]. His works established the mathematic base of modeling and algorithm on cable-structures.

Nonsmooth mechanics allows a general theoretical description but relevant problems must be solved numerically. In spite of many value contributions to the numerics of non-continuous systems [START_REF] Murty | Linear complementarity, linear and nonlinear programming[END_REF] the existing algorithms are still extremely time-consuming. Therefore the numerical solution of all kinds of complementarity problems is a topic of current research. Recently, many researchers are interested in the study of numerical algorithms subjected to unilateral constraints. The significant step may be to decide the real one of possible combination of constraints. For example, by means of common numerical method, if the number of cable members is equal to m, one must search for the real one among 2 m possible combination of constraints and therefore it consumes a large of time on computer. This is not suitable for calculating of a real construction. Also in this paper development of general computer method will be devoted to calculate nonlinear vibrations of cable-structures.

Actually, it has been noticed that non-smooth models overstep the framework of mechanical systems, since they also can apply for instance to electrical and NNW. In the present paper, the cable-structures are considered as a class of mechanical complementary-slackness systems. Based on the optimization algorithms for multi-body dynamics with unilateral contacts, an algorithm by means of artificial NNW has been developed. The following two classes of cable-structures have been considered: force-elongation of cable member follows elastic behavior and work-hardening assumption. Due to simplicity the former is used to prove the method reliability, and the latter as general cable-structure problem is handled. First, the complementarity problems for those structures have been formulated; then using generalized Gaussian' least action principle they are summarized as an optimization problem. Based on Hopfield's work, an artificial NNW has been designed and used to decide combination of possible constraints at each step in simulation. As examples, two cable-structures have been investigated. The calculated results for a simple suspension structure evidence the reliability and time-economization of the proposed method. An example of guyed mast shows the suitability of the proposed method for practical cable-structures.

Modeling of cable-structures as some complementarity problems

Although our work is aimed at cable-structures, the method proposed here is also available for some hybrid dynamical systems, like convex systems [START_REF] Ten Dam | Representations of dynamical systems described by behavioral inequalities[END_REF][START_REF] Ten Dam | The contact problem for linear continuous-time dynamical systems: a geometric approach[END_REF], or so-called complementary-slackness systems [START_REF] Van Der Schaft | The complementarity-slackness class of hybrid systems[END_REF][START_REF] Van Der Schaft | Complementarity modeling of hybrid systems[END_REF]. As examples in this paper, let us consider the following two systems with different nonlinear behaviors, to model those systems as complementary problems.

It is necessary to have a brief review of cable behaviors. It was concluded that the assumption of elastic behavior (Fig. 1a) of the structure yields the largest tensions in the cables [START_REF] Panagiotopoulos | A variational inequality approach to the inelastic stress-unilateral analysis of cable-structures[END_REF]. The elastic-prefectly plastic assumption (Fig. 1b) gives the lowest tensions in the cables. The work-hardening assumption (Fig. 1c) gives tensions between them. It is clear that the developed theory permits us to analyze the calble-structures in their "unstable region" because the possibility that a cable member can become slack is taken into account. In this paper, two cases, elastic behavior and work-hardening assumption, are considered.

The elastic assumption

The elastic-prefectly plastic assumption

Modeling of cable-structures as linear complementarity systems

It is well known that cable-structures can be modeled as linear complementarity systems [START_REF] Maier | Energy approach inelastic cable-structure analysis[END_REF][START_REF] Panagiotopoulos | A variational inequality approach to the inelastic stress-unilateral analysis of cable-structures[END_REF], but very little work for dynamical problems can be found. In the first example, a simple suspension structure is modeled as a Linear Complementarity Dynamical System.

Linear complementarity dynamical system

The governing equations of a Linear Complementarity Dynamical Systems [START_REF] Brogliato | Nonsmooth mechanics[END_REF][START_REF] Petersen | Chaotische Taumelschwingungen abgespannter Master[END_REF][START_REF] Pfeiffer | Multibody dynamics with unilateral contacts[END_REF] are written as

   ẋ = Ax + Bλ, y = C x + Dλ, y(t) 0, λ(t) 0, y(t) T λ(t) = 0.
(1) The modes of this hybrid dynamical system (there are 2 m modes, y and λ ∈ R m ) correspond to y i = 0, f or i ∈ I ⊆ {1, . . . , m}, λ i = 0, f or i ∈ I C , where I C is the complement of I in{1, . . . , m}. It is assumed that each mode is autonomous, which means that each consistent state x 0 there corresponds a unique couple(x(t), λ(t)), with λ(•)smooth.

If cable-elements are elastic, see Fig. 1a, some cable-structures can be described by Eqs. (1). To prove the reliability of the method proposed in this paper, a simple example has been considered firstly in the following subsection.

A simple example for a suspension structure

A simple suspension structure showing in Fig. 2 consists of a main cable, three slings and a beam. An unsymmetrical external force p(t)acts at l/4 position of the beam. Each sling is assumed to be elastic, see Fig. 1a.

Let the axial elongation of the i-th suspension sling be

l i = l i0 -w i -u i (i = 1, 2, 3), (2)
where l i0 -static expansion of the i-th sling; w i -elastic deflection of the beam at connected point with the i-th sling; u i -elastic deflection of the main cable at connected point with the i-th sling; The governing equations of the entire system are given as follows:

       m C ü -Hu = s i=1 F i δ(x -x i ), m B ∂ 2 w ∂t 2 + EI ∂ 4 w ∂ x 4 + CI ∂ 5 w ∂t∂ x 4 = s i=1 F i δ(x -x i ) + P(t), (3)
where m c and m B are mass permit unit length of the cable and the beam respectively. E is the Young's modulus of the beam; C indicates the damping factor; and P(t)describes the external excitation; F i is the force of the i-th sling acting at beam and main cable.

The boundary condition can be expressed by

u(0, t) = u(l, t) = 0; w(0, t) = w(l, t) = ∂ 2 w ∂ x 2 | x=0 = ∂ 2 w ∂ x 2 | x=l = 0.
Introducing the following relations

u = ∞ n=1 ψ 2n (x)v 2n (t); w = ∞ j=1 ϕ j (x)q j (t). ( 4
)
Take into account the boundary conditions, one has

φ j (x) = sin jπ l x (j = 1, 2, 3...); ψ 2n (x) = sin 2nπ l x (n = 1, 2, 3, ...). ( 5
)
Substituting expression (4) into Eq. ( 3), and integrating the resulting equations, the equations of motion for the entire structure can be rewritten as

       q j + 2ξ B j ω B j q j + ω 2 B j q j = 3 i=1 φ j (x i ) M B j F si + P j M B j , v2n + ω 2 C2n v 2n = 3 i=1 ψ 2n (x i ) M C2n F si , (6)
where

M B j = m B l 0 φ 2 j (x)dx, ω B j = jπ l 2 EI m B , ξ B j = Cω B j 2E , P j = l 0 φ j (x)Pdx; M C 2n = m C l 0 ψ 2 2n (x)dx, ω C 2n = 2nπ l H m C . Let Y = [Y q , Y v ] T = q 1 q 2 ...q j v 21 v 22 .
..v 2n T , the above equations can be rewritten in the following form:

Ÿ + Ẏ + Y = F K S + F P = F K K | L|λ + F P , S = K L = K (L -L 0 ) = K DY 0, (7)
in which,

L =   l 1 l 2 l 3   ; = B j 0 n ; = B j C n ; D = [ϕ j ψ j ]; F K =              ϕ 1 (x 1 ) M B 1 ϕ 1 (x 2 ) M B 1 ϕ 1 (x 3 ) M B 1 . . . . . . . . . ϕ j (x 1 ) M B j ϕ j (x 2 ) M B j ϕ j (x 3 ) M B j ψ 21 (x 1 ) M C 21 ψ 21 (x 2 ) M C 21 ψ 21 (x 3 ) M C 21 . . . . . . . . . ψ 2n (x 1 ) M C 2n ψ 2n (x 2 ) M C 2n ψ 2n (x 3 ) M C 2n              ; F P =       P 1 M B 1
. . .

P j M B j 0 n       ; K =   k 1 0 0 0 k 2 0 0 0 k 3   ; S = |S|λ; λ =    λ 1 λ 2 λ 3    ; |S| =   |S 1 | 0 0 0 |S 2 | 0 0 0 |S 3 |   ,
where λ i indicates as an indicative function of the i-th cable tension. According to the elastic tension-elongation relation, tension increment of each sling can be described in the following form.

S i = k i l i . ( 8
)
Let y i = s ik i l i and introducing the state vector,x = [Y Ẏ ] T , the mathematical problem (7) can be reformed into a standard Linear Complementarity Dynamical System as described in Eq. ( 2), in the following form.

ẋ = Ax + FP + FK S, y = S -K L S 0, y 0, S T y = 0 , ( 9
)
where

A = 0 I - - ; FP = 0 F P ; FK = 0 F K .
In Eqs. ( 9), one can see that i f S i 0, y i = 0,

I f S i = 0, l i < 0, y i > 0.
2.2 Modeling of cable-structures as a Mechanical Complementary-Slackness System

Formulation of Complementarity-Slackness Systems

If cables behave as the work-haderning, the cable-structures may be modeled as a class of Complementarity-Slackness Systems. Let us consider the second example, a guyed mast as a Complementary-Slackness System. This class of systems war defined in [START_REF] Van Der Schaft | The complementarity-slackness class of hybrid systems[END_REF][START_REF] Van Der Schaft | Complementarity modeling of hybrid systems[END_REF], namely

ẋ = f (x, u), y = h(x, u), y 0, u 0, y T u = 0. ( 10
)

Example 2: a guyed mast

Petersen had investigated a guyed mast [START_REF] Petersen | Chaotische Taumelschwingungen abgespannter Master[END_REF], which is described by DAE. In this mast the space chaotic vibration subjected to wind was discovered. But he did not consider the loosing phenomenon of cables, which may occur under the dynamic wind load. In this paper, based on Petersen's model a new model will be suggested as a Complementarity-Slackness System.

(1) Petersen's model The equations of space vibration of a guyed mast are given as follows [START_REF] Petersen | Chaotische Taumelschwingungen abgespannter Master[END_REF]:

m ẅx + d 1 ẇx + d 2 ẇ2 x + R x = F x (t), m ẅy + d 1 ẇy + d 2 ẇ2 y + R y = F y (t), ( 11
)
where m is effective mass, d 1 and d 2 are damping factors, R x and R y indicate the components of restoring forces in the x-and the y-directions; w x and w y describe the components of displacements at mast top in the x-and the y-directions; F x and F y are the components of wind load in the x-and the y-directions, respectively. Figure 3 shows the structure diagram of a simple guyed mast in which the x indicates the

(a) (b)

Fig. 3 The structure of a guyed mast direction of symmetric area on this mast, and r is distance from the mast to anchoring points. Introducing the dimensionless variables ξ and η, and let ξ = w x /r, and η = w y /r, the Eqs. ( 11) can be rewritten in the following dimensionless form:

ξ + d 1 m ξ + d 2 r m ξ 2 + S1 cos ᾱ -S2 cos β + S3 cos γ /rm = F x (t)/rm, η + d 1 m η + d 2 r m η2 + -S1 sin ᾱ + S2 sin β + S3 sin γ /rm = F y (t)/rm, ( 12
)
in which,

Si = S i cos δ, i = 1, 2, 3, ( 13
)
where Si is the horizontal projection of tension S i of the i-th cable, δ is defined as an inclined angle between the cable and the horizontal plan, and tension S i is the solution of the following algebraic equation:

S i S 0 3 + c i S i S 0 2 -a = 0, ( 14
)
in the above equation, the parameters are given by

a = g 2 s 2 cos 2 δ 24S 0 E s A s , b i = E s A s cos δr S 0 s , c i = a -b i i -1; i = i r , ( 15
)
where g is the weight per unit length of the cable; s means the length of the cable; E s stands for the Young's modulus; A s is cross-sectional area of the cable; S 0 describes pre-tensile force; i indicates the elongation of the i-th cable and its the first order approximations are given as follows:

1 = 1 + ξ - √ 3η + ξ 2 + η 2 -1 r, 2 = 1 -2ξ + ξ 2 + η 2 -1 r, 3 = 1 + ξ + 3η + ξ 2 + η 2 -1 r, ( 16
)
where i indicates dimensionless elongation of the i-th cable. In Eq. ( 12), ᾱ, β, γ are defined as inclined angles between the three cables and the x-direction, respectively, which can be determined using the following relations.

ᾱ = 60 • -arccos   2 + ξ - √ 3η 2 1 + ξ - √ 3η + ξ 2 + η 2   , β = arccos 1 -ξ 1 -2ξ + ξ 2 + η 2 , γ = 60 • -arccos   2 + ξ + √ 3η 2 1 + ξ + √ 3η + ξ 2 + η 2   .
(2) A new model as a Complementarity-Snackness System Some cables may loosen under a largest wind load. A new model is proposed here. It includes the nonsmooth behavior of the structure. Setting q = ξ η T and writing the Eq. ( 14) as S i = G( ), Eqs. [START_REF] Rockafellar | Convex analysis[END_REF] may be rewritten in the following form,

q + h( q) -Q(t) + S = q + h( q) -Q(t) + |S|λ = 0, S = G( ) -S 0 0, (17)
in which, λ i is the indicative function of the i-th cable, and where

h( q) = d 1 m q1 + d 2 r m q2 1 d 1 m q2 + d 2 r m q2 2 ; Q(t) = F x rm F y rm ; = cos ᾱ -cos β cos γ -sin ᾱ sin β sin γ ; S =   S1 0 0 0 S2 0 0 0 S3   ; λ =   λ 1 λ 2 λ 3   .
From the Eq. ( 14), the relationS = G( ) can be derived. Equation ( 17) can be rewritten as a standard form of a Complementarity-Snackness System described by Eq. [START_REF] Moreau | Application of convex analysis to some problems of dry friction[END_REF]. Let x = [q q] T and y = S -G( ), one has the following equations:

ẋ = f (x, S), y = S -G( (x)), y 0, S 0, y T S = 0. (18)
In the Eq. ( 18), it is clear that i f S 0, y = 0 and if S = 0, 0, y > 0.

Summing Eqs. ( 9) and ( 18), a general model for cable-structures may be expressed in the following form.

q + g(q, q) = H (q, q)S + P(t), y = Sh(q), S 0, y 0, S T y = 0.

According to Eq. ( 10), the above problem can belong to Compelementarity-Slackness Systems. It is well known that this problem can be solved using the optimization techniques. But some algorithms of optimization are still CPU-time consuming to decide possible constraints. Hence, a new algorithm of artificial NNW will be proposed in the following section.

The artificial NNW algorithm

Study of numerical algorithm to systems subjected to unilateral constraints is our main assignment. As a powerful tool by investigation of the dynamic systems with unilateral constraints, the optimization method has been commonly used to solve linear compelementarity problems (LCP) [START_REF] Pfeiffer | Multibody dynamics with unilateral contacts[END_REF]. Using generalized Gaussian' least action principle the Complementarity-Slackness Problem [START_REF] Petersen | Chaotische Taumelschwingungen abgespannter Master[END_REF] will be summarized as an optimization one.

It is generally agreed that the optimization problem can be solved by means of artificial NNW. The early work for analysis of LCP using the artificial NNW algorithm to solving LCP had been reported in [START_REF] Li-Zhi | A neural network for the linear complementarity problem[END_REF]. A new NNW based on a reformation of the LCP into unconstrained minimization problem was proposed. In this paper, based on Hopfield's work a class of artificial NNW will be designed and used to decide the combination of possible constraints at each step of simulation.

According to Hopfield's method, we put ∂ E/∂ V i = -dx i /dt. Whent → ∞, and dx i /dt = 0, a stable solution can be obtained as

   X j+1 = W j V j + I j , V j+1 = F X j = 0, X j 0 1, X j > 0 , ( 23
)
in which, W j = -∇ q f (q j ) T H (q j , q j )|S|,

I j = -[∇ q f (q j ) T z j + d dt (∇ q f (q j )) T q j ]
= -[∇ q f (q j ) T (-g(q j , q j ) + P j ) + d dt (∇ q f (q j )) T q j ].

Using artificial NNW to decide the real involved constraints and by means of a modified Rung-Kutta method, the solution of Eqs. ( 21) can be obtained. The flow chart for this method is shown in Fig. 4. The description of the flow chart:

(1) Using Runge-kutta method, response of the whole structure given by (19, a) at the j-1-step can be obtained; (2) Introducing the above results into equation of artificial NNW [START_REF] Loetstedt | Mechanical systems of rigid bodies subject to unilateral constraints[END_REF], the values λ at the j-th step can be determined and the involved constraint combination at the next step can be decided; (3) Substituting the values λ at the j-th step into the Eq. (19, a) and return to Eq. ( 23).

Numerical examples and discussion

Example 1: a simple suspension structure

Figure 2 shows a simple suspension structure, which consists of a main cable, three slings and a beam. The structure parameters are assumed as listed below: length of beam l = 8 m, bending stiffness EI = 1.2 × 10 7 N/m 2 , mass per unit length of beam m b = 650 kg/m, stiffness of the sling k = 3.75 × 10 3 N/m, the three slings are arrange at l/4, l/2, 3l/4 of beam; mass per unit length of the main cable m c = 100 kg/m, the droop f = 0.1 m; the external excitationP = P 0 sin(θ t) with P 0 = 5, 000 N, θ = 0.75ω B 1 , acting at l/4 of the beam.

Figures 5 and6 exhibit the time histories of the deformations and tensions for the three slings using both methods of artificial NNW and general search method (GSM), respectively. In those figures, the results calculated by the both methods are the same. This testifies that the method of artificial NNW is reliable. Moreover, the calculated results show that using the NNW method to solve this simple problem requires only 15 min of CPU-time, but using GSM needs several hours. From the above mentions one can conclude that it may be better, to use the artificial NNW by solving the dynamics of the system with unilateral constraints. Figures 7 and Fig. 8 give the results by using Petersen's model and by new model suggested in this paper, respectively. In Fig. 7, one can find negative values of tension, this is not reasonable, while this phenomenon is avoided in the modified model. In Fig. 8, some cables sometime may be loosening, and their tensions are equal to zero. It may be interpreted, that for this example the loosening phenomenon occurs due to the insufficient pre-tensions in cables. From the calculated results, it can be discovered that the three cables can be loosening, but they can be loosening not at the same time. In other words, no restoring forces in the motion equations are vanished, so that the mast does not collapse. However, because of non-smoothness, this mast exhibits complex dynamical behavior. Hence, this mast may be damaged due to fatigue.

Conclusions

Through the analysis of two examples, the following conclusions may be drawn:

(1) The loosening phenomenon of some cables may occur when the design of structures with cables is not reasonable. (2) In modeling of cable-structures, consideration of the non-smooth behavior of those systems has been found to be necessary. [START_REF] Murray | Analysis of ineladtic suspension structures[END_REF] The artificial NNWs used to analyze the dynamic systems with unilateral constraints may be a reliable and suitable method, which can save a large of time on computer.

Fig. 1

 1 Fig. 1 The elongation-tension curves of cables. a The elastic assumption. b The elastic-prefectly assumption. c The workhardening assumption

Fig. 2 A

 2 Fig. 2 A simple suspension structure

4. 2 Fig. 4

 24 Fig.[START_REF] Tene | Nonlinear space trusses with tension-only members[END_REF] The flow chart for the method of artificial neural network

Fig. 5

 5 Fig. 5 The time histories of the tensions of the slings

Fig. 6

 6 Fig. 6 The time histories of the deformations of the slings

Fig. 7

 7 Fig. 7 The results calculated using Petersen's model. a The tension of cable 1 calculated using Petersen's model. b The tension of cable 2 calculated using Petersen's model. c The tension of cable 3 calculated using Petersen's model

Fig. 8

 8 Fig. 8 The calculated results using new model. a The tension of cable 1 calculated using new model. b The tension of cable 2 calculated using new model. c The tension of cable 3 calculated using new model

The generalized Gaussian's least action principle

The Gauss's least action principle is one of the several variational principles in mechanics. In fact, Gauss's principle is a reinterpretation of d'Alember's principle into a minimum principle. It turns out that Gauss's principle is also valid for systems subjected to unilateral constraints [START_REF] Bolotov | On Gauss' principle[END_REF][START_REF] Loetstedt | Mechanical systems of rigid bodies subject to unilateral constraints[END_REF][START_REF] Moreau | Quadratic programming in mechanics: dynamics of one sided constraints[END_REF][START_REF] Sinitsyn | The principle of least constraint for systems with non-restoring constraints[END_REF]. It is found that integrating of the system at a given time instant is a quadratic procedure, and is closely related to complementarity conditions. The general form of the complementarity conditions can be found in the works by Loetsted [START_REF] Loetstedt | Mechanical systems of rigid bodies subject to unilateral constraints[END_REF], Baraff [26] and Moreau [START_REF] Moreau | Quadratic programming in mechanics: dynamics of one sided constraints[END_REF]. From the above dynamics one may see that a solution to the LCP ( 19) is equivalent to the following quadratic programming problem [START_REF] Pfeiffer | Multibody dynamics with unilateral contacts[END_REF], and vice-versa.

z = -g(q, q) + P(t),

From this, the following dynamical equation can be obtained.

q + g(q, q) = H (q, q)S + P(t) = H (q, q)|S|λ + P(t), S = f (q) 0, [START_REF] Li-Zhi | A neural network for the linear complementarity problem[END_REF] where the inequalities S = f (q) 0 can be rewritten into the following complementarity conditions: y = Sf (q) 0; S 0; S T y = 0.

Thus it is proved that the optimization problem (20) and the dynamics with unilateral constraints [START_REF] Petersen | Chaotische Taumelschwingungen abgespannter Master[END_REF] are equivalent to each other.

An algorithm of artificial NNWs

Much of the resurgence of interest in NNWs during the early 1980s can be attributed to the work of John Hopfield [START_REF] Hagan | Neural network design[END_REF]. He wrote two highly influential papers in 1982 [START_REF] Hopfield | Neural networks and physical system with emergent collective computational properties[END_REF] and 1984 [START_REF] Hopfield | Neurons with graded response have collective computational properties like those of two-state neurons[END_REF], respectively. A number of important ideas and a clear mathematical analysis in his early papers include content-addressable memory, analog-to-digital conversion, and optimization [START_REF] Hopfield | Neural' computation of decisions in optimization problems[END_REF][START_REF] Tank | Simple 'neural' optimization networks: an A/D converter, signal decision circuit and a linear programming circuit[END_REF]. The Hopfield network does not have a learning law associated with it. It is not trained, nor does it learn on its own. Instead, a design procedure based on the Lyapunov function is used to determine the weight matrix and the bias vector. This trick, of course, is of convenience.

In this section, adopting Hopfield's method [START_REF] Hopfield | Neural' computation of decisions in optimization problems[END_REF], an artificial NNWs equivalent to optimization problem (20) will be designed. Let V i = λ i , V i = F(x i ), from the Gaussian's function, an energy function may be assumed to be

V T ∇ q f (q) T H (q, q)|S|V + V T (∇ q f (q) T z )