Jürgen Schnerr
email: jschnerr@fzi.de

Oliver Bringmann
email: bringmann@fzi.de

Wolfgang Rosenstiel
email: rosenstiel@fzi.de

Cycle Accurate Binary Translation for Simulation Acceleration in Rapid Prototyping of SoCs *

In this paper, the application of a cycle accurate binary translator for rapid prototyping of SoCs will be presented. This translator generates code to run on a rapid prototyping system consisting of a VLIW processor and FPGAs. The generated code is annotated with information that triggers cycle generation for the hardware in parallel to the execution of the translated program. The VLIW processor executes the translated program whereas the FPGAs contain the hardware for the parallel cycle generation and the bus interface that adapts the bus of the VLIW processor to the SoC bus of the emulated processor core.

Introduction

Nowadays, semiconductor chips reach such high complexities, that the integration of complete systems with high functionality on a single chip (so called Systems-on-Chip, SoC) becomes possible. Using an SoC bus, these SoC components are connected on the chip.

The development of software for such SoCs is problematic as the software running on SoCs usually includes very hardware near parts like device drivers. The requirements for this software are that I/O accesses to the bus must be cycle accurate in order to make it possible to validate the bus interfaces to the hardware or the handshakes on the bus.

Existing solutions for this problem are the hardware/software cosimulation [START_REF] Rowson | Hardware/Software Co-Simulation[END_REF] of processors and the attached hardware on the bus or the emulation of the processor using a hardware emulator. Using hardware/software co-simulation, the most accurate way of simulation would be a simulation of the processor and the attached hardware on an HDL simulator. But this leads to very slow execution times. To speed up the simulation, an abstraction of the processor has to take place. This can usually be done in two different ways. Either the source code has to be recompiled for the processor of the machine the simulation is running on, or an instruction set simulation (ISS) has to be used. An ISS can be implemented as an interpreted ISS, a just-in-time compiled ISS, or as a compiled ISS. Using a bus functional model (BFM), the recompiled code or the ISS has to be coupled to an HDL simulator simulating the hardware attached to the processor. These existing solutions using hardware/software co-simulation suffer from different disadvantages, especially concerning accuracy and/or execution speed. Further disadvantages are, that some solutions require the source code of the program. Also if the hardware has to be simulated, a model of the hardware environment is needed.

Using hardware emulators or FPGAs for the emulation of the processor can lead to a significant speed-up of execution speed [START_REF] Schmitt | Verification of a Microcontroller IP for System-on-a-Chip Designs Using Low-Cost Prototyping Environments[END_REF]. The main advantage of such solutions is that a very accurate emulation of the processor is possible at a reasonable speed. The disadvantage is that the detail level of the emulation cannot be chosen, therefore it can be more detailed than needed. Also a synthesizable specification of the processor is needed.

The rapid prototyping system used in this paper was presented in [START_REF] Schnerr | Instruction Set Emulation for Rapid Prototyping of SoCs[END_REF]. It is based on the hardware presented in [START_REF] Haug | VLIW Based Emulation of Digital Designs with the RAVE System[END_REF][START_REF] Haug | A Hardware Platform for VLIW Based Emulation of Digital Designs[END_REF]. This system was designed to overcome the disadvantages of previous existing systems by offering a fast execution speed of the simulated system with sufficient cycle accuracy. To find a trade-off between execution speed and accuracy, several detail levels of code execution are offered. Another advantage compared to other solutions is, that neither a HDL description of the processor nor source code of the running program is needed.

Basically, the presented prototyping system consists of two parts. One part is the TMS320C6201 (C6x) VLIW processor from Texas Instruments [START_REF]TMS320C62xx CPU and Instruction Set Reference Guide[END_REF], the other part consists of FPGAs. The VLIW processor is used for the emulation of the instruction set of the SoC processor core. As examples of processor cores the TriCore processor from Infineon [START_REF]TriCore Architecture v1.3 Manual[END_REF] and the ARM processor [START_REF]ARM7 Data Sheet[END_REF] were used. This paper will focus on the emulation of the processor core, that we have implemented using a compiled ISS. The binary translator that carries out the compiled ISS will be described in more detail in this paper. Compared to existing binary translators, this one annotates the translated program with instructions that access the cycle generation hardware in the FPGAs. The FPGAs of the system are also used for a bus interface between the C6x processor and the attached hardware that expects to be connected to an SoC bus.

Basically, the cycle accurate translation works as follows: The runtime of each basic block of the untranslated program is predicted. At the beginning of the translated basic block, code is included that triggers the start of the cycle generation with the number of predicted cycles. When such a translated basic block is run, it starts the cycle generation with its first instruction, and then the rest of the basic block runs in parallel to the cycle generation. At the end of the basic block there is an instruction that synchronizes the cycle generation with the instruction execution before the execution of the next basic block begins.

This method of letting the translated program trigger the clock cycles for the hardware itself reduces the bottleneck between the processor executing the program and the attached hardware.

The remainder of this paper is structured as follows: Section 2 shows a short overview on the state of art in using recompiled source code and ISS under consideration of timing issues. Section 3 of this paper describes the implementation of the cycle accurate binary translator. Section 4 provides results. Section 5 concludes the paper and gives an outlook.

State of the art

As stated in the previous section, a recompilation of the source code or an ISS can be coupled to a BFM to provide a cycle accurate abstraction of a processor.

Using a recompilation of the source for the target processor usually changes the timing of the program completely. To include the timing of the original processor there are several possibilities. One possibility is to add timing information in the form of delay annotation to the source code [START_REF] Bacivarov | Timed HW-SW Cosimulation Using Native Execution of OS and Application SW[END_REF][START_REF] Yoo | Automatic Generation of Fast Timed Simulation Models for Operating Systems in SoC Design[END_REF]. Another approach to the inclusion of timing information was shown in [START_REF] Lee | Timed Compiled-Code Simulation of Embedded Software for Performance Analysis of SOC Design[END_REF]. In this case, the authors modified a compiler to analyze the program being translated and to include the timing information into the generated code.

When a recompilation of the source code is not possible or too inaccurate, an ISS can be used. There are three possibilities for the implementation of such an ISS. The first one is the interpretative simulation. At the moment, this is the most commonly used method. As the interpretative simulation has to decode the instructions during the runtime, it suffers from low performance. To overcome this problem, a just-in-time (JIT) compilation (also called dynamic compilation) of the program can be used [START_REF] Nohl | A Universal Technique for Fast and Flexible Instruction-Set Architecture Simulation[END_REF]. Here the performance is much better as the program is translated during the runtime and blocks of the translated program are cached leading to a much better performance if the same part of the program has to be executed more than once.

The third possibility is the compiled simulation (also called binary translation or static compilation) of the program. As the decoding and the translation of the program are done statically during the compiletime it can reach the fastest execution speed of the three possible implementations.

The development of binary translators started at the end of the 1980s where they were needed for the translation of existing CISC programs in order to migrate to RISC architectures. For this reason they offered only a functional translation of the program with no consideration of timing issues. A good introduction and an overview of such translators are shown in [START_REF] Cifuentes | Binary Translation: Static, Dynamic, Retargetable[END_REF].

For the translation of programs that access hardware these purely functional solutions were not enough. Therefore, Cogswell and Segall presented in [START_REF] Cogswell | Timing Insensitive Binary to Binary Translation[END_REF] a methodology of translating the program where the translated program is instrumented with timing information in a way that during the execution of the program on the target platform the timing for the source platform contained in the timing information and the timing on the target platform can be compared at regular intervals and the program execution can be delayed if it is necessary. This paper will introduce a cycle accurate binary translation specially tailored to the needs of prototyping of SoCs. In contrary to previous solutions, it works on a very low, hardware near level.

Cycle accurate static compilation

This section will describe the compiler presented in this paper in more detail. The design of the compiler is made up of two parts. One part consists of the description of the instruction set and the architecture of the processor whose code has to be translated (the source processor). This description is done using C++ classes. The other part is independent of the source processor and is implemented as a library. To generate a working compiler both parts have to be linked together. This way, the compiler can be adapted to different source processors.

Although it is possible to write the classes containing the description of the source processor by hand, it is not recommended due to its errorproneness. Therefore, this processor is usually defined in an XML file that is translated into the appropriate C++ code by a tool. This XML file contains an architecture description and a description of the instruction set of the processor. The architecture description contains a description of the pipelines and the caches of the processor, whereas the instruction set description contains information about instruction decoding and the semantics of the described instruction written in an intermediate code which resembles the assembler instructions of the C6x processor but does not have their constraints.

Figure 1 shows an overview of how the compiler works. The parts of the compiler needed for cycle accurate translation of code are shown in grey.

Figure 1. Cycle accurate static compiler

Using an appropriate class, the compiler reads the object file, which is usually provided in ELF format (other formats are also possible). Afterwards, this object code will be decoded and translated into an intermediate representation. This intermediate representation consists of a list of objects. Each of these objects is representing one intermediate instruction.

In a next step, the basic blocks of this program are found out using the list containing the translated program, and a list of basic blocks is built.

Using the list of basic blocks, the base addresses of load/store instructions have to be found out, as far as this is statically possible. This is needed for two reasons. One reason is to change the base addresses of load/store instructions accessing memory to the new memory addresses of the target system. The other reason is to find out, which of these load/store instructions are I/O instructions. These I/O instructions have to be replaced by instructions accessing the hardware of the bus model.

After that, a static calculation of the number of cycles each basic block would have taken on the source processor is made. How this calculation works is described in more detail in Section 3.3.

In the next step, the code for the cycle generation for each basic block has to be inserted into the basic blocks. How this code looks like and how it works is described in Section 3.1.

Not every effect of the processor architecture on the number of cycles can be predicted statically. Therefore, if effects of branch prediction and caches should be taken into account, additional code has to be added in this step. Further details concerning this code are described in Section 3.4.

After that, further transformations of the intermediate code have to be done. Here, instructions that can be executed in parallel on the VLIW processor have to be found out. Furthermore, software pipelining of the instructions and register binding have to take place. Also, every instruction has to be assigned to the functional unit it will run on.

The result is a program for the TMS320C6x VLIW processor annotated with cycle information. More information about how this annotation works is described in the next section.

Annotation of translated code

Usually, there is an external common clock for the processor and the attached hardware, but in this case there is no such common clock. Instead, the emulated processor has to generate the clock cycles for the attached hardware itself. This is done using accesses to a hardware device in the FPGAs called synchronization device. As shown in Figure 2 the compiler adds an instruction that starts the cycle generation at the beginning of the basic block. This instruction is a write access to the synchronization device that contains the number n of cycles this basic block would need on the source processor. From now on the execution of the instructions in the translated basic block and the generation of the cycles for the attached hardware run in parallel until the executed program reaches the "wait for end of cycle generation" instruction. This instruction consists of a read access to the synchronization device. If the cycle generation is already finished, this instruction does nothing. If it is not finished, this instruction waits until the generation is finished.

This way, the execution of the translated program and the attached hardware can be synchronized without having the bottleneck of permanent hardware accesses by the executed program.

Cycle calculation

In order to guarantee an as fast as possible execution of the translated code, static cycle calculation, as mentioned in the previous section, should be used if it is possible. But compared to a dynamic cycle calculation this leads to less accuracy as not every effect of processor architecture on the number of executed cycles can be exactly predicted.

In order to guarantee both an as fast as possible execution of the code as well as the highest possible accuracy, it is possible to choose different levels of accuracy of the generated code when the compiler is run.

There are the following three detail levels of cycle accuracy:

1. purely static prediction 2. dynamic improvement of the static prediction (modeling of the branch prediction)

dynamic inclusion of instruction caches (additionally to point 2)

The cycle calculation in these different detail levels will be discussed in more detail in the following text.

Static cycle calculation of a basic block

The simplest possible way of calculating the execution time a basic block consists of the summation of the execution or the latency times of the instructions of this basic block. But this is only possible with simple architectures. Concerning modern architectures, it is too inaccurate as it does not consider pipeline effects, super scalarity and caches of these architectures.

In order to predict pipeline effects and the effects of super scalarity statically, modeling the pipeline per basic block becomes necessary [START_REF] Lim | An Accurate Worst Case Timing Analysis for RISC Processors[END_REF]. Using information about the instructions and the pipelines of the processors, it can be found out which instructions of the basic block can be executed in parallel with a super scalar processor and which combinations of instructions in the basic block cause pipeline stalls.

With the information gained by the modeling of the basic block a prediction can be made how many cycles the basic block would have needed on the source processor.

In the following section, it will be shown how such a prediction can be improved during runtime.

Dynamic correction of cycle prediction

Not every cycle that a processor executes can be correctly predicted in a static way. For example, if a conditional branch occurs at the end of a basic block or the effects of instruction caches during the execution of a basic block.

In these cases, additional code that counts correction cycles has to be inserted into the translated basic block, and at the end of the basic block a correction block has to be inserted which generates these additional correction cycles before the next basic block can be run.

Figure 3 shows such an annotated basic block that offers additional dynamic correction of the cycle count. The additional cycle counting code for the consideration of branch prediction is explained in the next section. The division of the basic block for the calculation of additional cycles for instruction cache misses as shown in the Figure 3 wrong and whether the branch was taken or was not taken. Usually, such a conditional branch needs a minimum number of cycles in all cases. This number can be added to the cycle generation code of the basic block. Whether there are additional cycles needed or not, has to be found out using additional code which has to be inserted in this step just before the conditional branch. During runtime, the inserted code has to check whether the branch prediction was correct or not and whether the branch is taken or not. The additional cycles will be added to a cycle counter for the execution in the correction block.

Instruction cache

For additional simulation of the instruction cache, three steps have to be done. First, space has to be reserved at the end of the translated program. This space is used to hold the data of the simulated cache during runtime. The second step that has to be done, is that every basic block of the translated program has to be divided into cache analysis blocks. In a third step, code for the cache handling has to be added to each of these cache analysis blocks.

All three steps will be explained in detail in the following three sections.

Saving cache data At the end of the translated program space for cache data is added. This space holds the valid bit, the cache tag and the least recently used (lru) information (containing the replacement strategy) for each cache set during the runtime. To simplify the handling of the cache tag and the valid bit, they are combined into one word which is usually 32 bits long.

The number of cache tags and the according valid bits is depending on the associativity of the cache, e.g. for a two way set associative cache there are two of them.

Cache analysis blocks For the consideration of the effects of instruction caches, each basic block has to be divided into cache analysis blocks. Each cache analysis block contains that part of a basic block that fits into a single cache block. Every instruction of such a cache analysis block has the same tag and the same cache index.

To find out the cache analysis blocks of a basic block all addresses of the untranslated instructions of this basic block have to be examined. Each address consists of tag information, cache index and block offset information.

The cache index information determines at what position of the cache the instruction with this address is cached. The tag information has to be used to find out which address had been cached, as there can be multiple addresses with the same cache index. The block offset is the position in the cache line at which the instruction is cached. For this simulation, the block offset information is not needed, as no real caching of data takes place. Now that the cache analysis blocks have been found out, code that handles the simulated cache and calculates the additional cycles of cache misses has to be added to these blocks. How this is done is described in the next section.

Cycle calculation code During runtime, it has to be found out, for each cache analysis block the basic block consists of, whether it is in the simulated cache or not. This way, cache misses can be detected.

Each cache analysis block is characterized by a tag and cache set index combination. At the beginning of each cache analysis block there is a branch to a subroutine included. This subroutine, shown in Figure 4, gets the tag and cache set combination as parameters. It checks whether the tag of the cache analysis block can be found in the specified set and whether the valid bit for the found tag is set.

input: tag and cache set index of cache analysis block output: changed cache information and cycle correction counter if tag can be found in specified set and valid bit is set then renew least recently used information // cache hit else // cache miss use lru information to find out tag to overwrite write new tag set valid bit of written tag renew lru information add additional cycles to cycle correction counter return to normal execution of the cache analysis block

Figure 4. Cycle correction for caches

If yes, the block is already cached und no additional cycles are needed. Only the lru information has to be renewed.

If no, the lru information has to be used to find out which tag has to be overwritten. After that, the new tag has to be written instead of the found old one, and the valid bit for this tag has to be set. Also the lru information has to be renewed. In a last step, the additional cycles have to be added to the cycle correction counter.

After that, the subroutine returns, and the normal execution of the instructions of the cache analysis block can continue.

The code of the subroutine is automatically generated by the translator on basis of the description it has about the cache of the source processor. It is appended to the translated code. In large basic blocks, this code can be included into the basic block making the subroutine call unnecessary and the parallel execution of the cache calculation code and the executed program on the VLIW processor possible.

Debugging of translated code

The debugging of code annotated with cycle information is implemented using an interface program between the translated code and the remote debugging interface of the GNU Debugger (gdb). One important problem which had to be solved was the implementation of a single step possibility for the debugger as the translated code only contains cycle generation for complete basic blocks.

This has been solved in the way that the debug code contains two translations of the original code. In one of these translations the code has to be annotated with a basic block oriented cycle generation, and in the other one it has to be annotated with an instruction oriented cycle generation. The code with the instruction oriented generation of cycles contains the cycle generation for each translated instruction and a branch after each such instruction into the ROM routines that communicate with the debug interface.

Using the described two translations of the code, the debug interface has to implement break points, single step execution and normal program execution. Also using information made during the translation of the program, the debug interface has to translate the register names and the addresses used.

Break points that occur in a certain basic block are always set at the beginning of a basic block of the basic block oriented translated code. To get to the real break point the single step program has to be used.

Results

In order to test the execution speed and the accuracy of the translated code, a few examples were compiled using a C compiler into TriCore object code. This object code was then translated into cycle accurate TMS320C6x binary code and executed on the emulation platform. As a reference, the execution speed and the cycle count of the TriCore code has been measured on a TriCore TC10GP evaluation board, using the methodology described in [START_REF] Zheng | An Easy Way to Measure TriCore Cycle Performance[END_REF].

The examples consist of two more control flow dominated programs (gcd, sieve), two filters (fir, ellip), and two programs that are part of audio decoding routines (dpcm, subband). A fast execution speed of the cycle accurate code can be gained especially for examples with large basic blocks like ellip and subband. This is due to the fact that a program containing larger basic blocks needs fewer cycle generation instructions and fewer branches. Branches are very costly in regard to the cycle time on the C6x processor. Also larger basic blocks allow a better parallelization of the instructions they contain on the VLIW processor, leading to a smaller execution time.

On the other side, the program sieve contains many small basic blocks. Each one of these basic blocks has its own cycle generation code. Therefore, the translated code containing the cycle information is much slower than the translated code without information as can be seen in Figure 5. For the compiled code without cycle information the C6x processor needs about three cycles to execute a translated TriCore instruction. This relatively low overhead results from the fact that in the translated code of these examples on the average about two or three C6x instructions can be executed in parallel.

The addition of cycle information to the code adds on the average only a little more than one clock cycle per translated TriCore instruction. Also the overhead for the additional TriCore branch prediction is very low.

On the next detail level, the additional consideration of instruction caches, about six times more cycles are needed in comparison to the code containing the branch prediction. This is due to the fact that for each cache analysis block a complete routine has to be run to find out whether there is a cache miss or a cache hit. Table 2 provides a comparison with the results published in [START_REF] Schmitt | Verification of a Microcontroller IP for System-on-a-Chip Designs Using Low-Cost Prototyping Environments[END_REF]. It uses the same three example programs and compares the translation results with the results running the programs on an RT level simulation of the TriCore processor core on a workstation and the result of the emulation of the core on a prototyping platform using a Xilinx XCV2000E FPGA. The emulation on this FPGA is running with a frequency of 8 MHz.

The table shows that the translation results considering caches are about in the same range as the emulation results on the FPGA. Whereas the translation results of the other two detail levels are significantly (28 and 42 times) faster for the sieve example and about 3 to 5 times faster for gcd and fibonacci. The relatively slow execution of the latter two examples results from the very short basic blocks they consist of.

Apart from the possibility to gain a much higher performance by lowering the detail level, the translation approach has another important advantage compared to the emulation of the processor core with an FPGA. As processors get more and more complex and as there is a capacity gap between the ASIC and FPGA gate capacities, it will remain problematic to map state-of-the-art processors onto FPGAs. Using the translation approach, a more complex processor only results into a more complex architecture and instruction set description for the translator.

Conclusion and outlook

In this paper, we have shown how a microprocessor can be emulated cycle accurately on a prototyping platform using the annotation of translated code.

Especially pipelines and caches had to be considered. It has been shown, how cycle generation code for these can be generated and inserted. This inserted code combined with the appropriate hardware adds only a low overhead to the translated code and allows a fast execution. To find a trade-off between accuracy and execution speed, the cycle accuracy of the translated code can have several detail levels. In addition to the possibility to choose the detail level of accuracy, another important advantage of the described emulation system is that no RT level code of the emulated processor has to be available.

There are several possibilities to enhance the system to make it either faster or more accurate. Execution speed improvements can be achieved by the addition of the possibility of cycle generation for blocks which are larger than the basic blocks that are used now. Another speed improvement can be, to combine the translation of the code with a static cache analysis. This could speed up the translated code regarding the simulation of caches. To improve the accuracy of the translated code, additional consideration of data caches can be added. Furthermore, the accuracy of certain instructions could be improved. A consideration of data dependent behavior of certain instructions like multiplications and divisions can also help to improve the accuracy in certain cases. For example, on a processor that uses a Booth multiplier the delay of this multiplier depends on operand value. Also, memory access delay should be taken into account for load/store instructions. Such a delay can happen if, for instance, a cache miss happens and both the processor (via cache) and a DMA controller want to access the memory.

Another enhancement of the system that is especially important for the execution of embedded software and the execution of operating systems will be the consideration of interrupt handling and exceptions.

  program further transformations of intermediate code object file insertion of cycle generation code insertion of dynamic correction code

Figure 2 .

 2 Figure 2. Annotated basic block

Figure 3 . Annotated basic block with dynamic correction 3 . 4 . 1 .

 3341 Figure 3. Annotated basic block with dynamic correction 3.4.1.Branch prediction Conditional branches usually have different cycle times considering whether the branch prediction was right or wrong and whether the branch was taken or was not taken. Usually, such a conditional branch needs a minimum number of cycles in all cases. This number can be added to the cycle generation code of the basic block.Whether there are additional cycles needed or not, has to be found out using additional code which has to be inserted in this step just before the conditional branch. During runtime, the inserted code has to check whether the branch prediction was correct or not and whether the branch is taken or not. The additional cycles will be added to a cycle counter for the execution in the correction block.

Figure 5 .

 5 Figure 5. Comparison of speed

Figure 6 .

 6 Figure 6. Comparison of cycle accuracy

Table 1 . Clock cycles per TriCore instruction

 1 The average number of clock cycles needed for the execution of one TriCore instruction is shown in Table1. The results on the C6x processor are the average value of all examples.

	TC10GP Evaluation Board	1.08
	C6x without cycle information	2.94
	C6x with cycle information	4.28
	C6x branch prediction	5.87
	C6x caches	35.34

Table 2 .

 2 3 % for the program Software runtime comparison ellip to 15 % for the program sieve. Especially for control flow oriented programs like gcd it can be seen that branch prediction can play an important role.

		gcd	fibonacci	sieve
	# of executed instructions	1484	41419	20779
	Simulation	approx.	approx.	approx.
	(Workstation)	28 sec	10 min	18 min
	Emulation (FPGA)	321 µsec 3.9 msec 21.8 msec
	Translation
	C6x cycle	63.1 µsec 950 µsec	520 µsec
	C6x branch	94.6 µsec 1.4 msec	781 µsec
	C6x cache	416 µsec 6.3 msec	5 msec

Proceedings of the Design, Automation and Test in Europe Conference and Exhibition (DATE'05) 1530-1591/05 $ 20.00 IEEE

Proceedings of the Design, Automation and Test in Europe Conference and Exhibition (DATE'05)

* This work has been partially supported by the BMBF/MEDEA+ project SpeAC A508 under grant number 01M 3067 A.


Nicolas Holvoet
email: nicolas.holvoet@alumni.enac.fr

Wasanchai Vongsantivanich
email: wasanchaiv@gistda.or.th

Supatcha Chaimatanan
email: supatcha@gistda.or.th

Daniel Delahaye
email: delahaye@recherche.enac.fr

Mission planning for a non-homogeneous Earth observation satellite constellation for disaster response

This article explores the mission planning problem for a constellation of non-homogeneous Earth observation satellites for disaster response to promote cooperation in the Asia-Pacific region between different space agencies or satellite operators. Given a set of weighted requests consisting of a polygonal area, a list of collaborating satellites and a time horizon, the goal is first to compute the observation opportunities and then to plan a sequence of acquisitions for each sensor to optimize the total coverage in the first days to help the disaster management agencies. The non-homogeneous satellites (non-agile or agile) and the different orbits result in a huge combinatorial search space. A simplified computational method is proposed to generate the observation opportunities and a metaheuristic based on the Simulated Annealing algorithm is used to solve the mission planning problem. Real requests from disaster management agencies are used to evaluate the proposed methodology.

I. Introduction

S ace-based Earth observation has become a common practice to get invaluable information about Earth and human activities. For disaster mitigation and disaster management, prompt remote sensing data from space covering the affected area is indispensable for planning, analyzing the situation for disaster relief, post-disaster management or even future disaster prevention. Having more data, both spatially and temporally, will undeniably lead to clearer situation assessment and more effective decision-making.

Considering using a single remote sensing satellite, there are significant limitations in spatial coverage, temporal resolution and varieties of data to assimilate to perform further analysis. With the current global trends in space utilization, it can be seen that the number of satellites launched into orbit has increased progressively every year. To maximize the benefits from these space resources, constellations of satellites are designed. However, most of the Earth observation satellites (EOS) that are currently working in space are individually operated by different countries or agencies. Therefore, acquisition plan for each satellite is also individually performed based on each acquisition request that is submitted to each particular operator. For this reason, there are many initiatives to collaborate and coordinate satellite utilization among the satellite operators, for example: CEOS's Virtual Constellation, International Charter, Disaster Monitoring Constellation, etc. "Sentinel Asia" is one of these initiatives, established in 2005 to promote international cooperation to prevent, mitigate and respond to natural disaster damage in the Asia-Pacific region by using space technology. In case of a major event, it aims to allow collaboration between space agencies to form a constellation of different types of Earth observation satellites which will perform a mission together.

The Geo-Informatics and Space Technology Development Agency (GISTDA) has developed a platform called OPTEMIS (Operation Planning Tool for Earth observation MISsion) to support the Sentinel Asia initiative. The goal of this tool is to facilitate and improve the operational workflow of the space Data Provider Nodes (DPN). This platform will be used as an intermediary to collect imaging requests from different users, perform feasibility assessment, manage the mission plan and track the mission progress, which will help in synergizing closer collaboration between the stakeholders (rescue and disaster response teams, Data Provider Nodes, Data Analysis Nodes, program coordinator, etc.).

We present below some general considerations concerning the problem, including the goals and the main constraints. Section II gives an overview of the previous work done on similar problems. The model and the acquisition feasibility assessment module are introduced in section III. We present in section IV the mission planning and task allocation algorithm that is used to compute acquisition requests for each collaborating satellite. Finally, section V demonstrates the benefit of the proposed system on real-world scenarios.

A. Agile satellites

In this paper, we consider a constellation of non-homogeneous satellites. Non-agile satellites have no maneuvering ability and the area they can acquire is determined by their instantaneous position on the orbital trajectory. Semi-agile satellites have one degree of freedom on the roll axis: the instrument can rotate to extend its field of regard (FOR) on the left and/or on the right side of the satellite's trajectory. This oblique viewing ability increases the imaging opportunities of a given area during an orbital cycle by reducing its revisit delay. Agile satellites enable off-nadir access by moving along their roll axis like semi-agile satellites but also along their pitch axis to point ahead or backward. Thus, the starting date of the acquisition can be freely decided in a time interval which depends on the minimum and maximum pitch angles allowed. This extended time window allows to schedule more observations which would not always have been feasible during a given orbit by a non-agile or semi-agile satellite as shown in Figure 1.

For both non-agile and (semi-)agile satellites, a minimum transition duration must be computed between two candidate observations to check if they can be successively acquired with enough time for the processing (compression, storage, downlink. . . ) and maneuvering actions (changes in roll and pitch angles). Naturally, image quality strongly depends on the pitch and roll angles and could be a hard or soft constraint from the user.

B. Problem statement

In case of a major disaster, authorized users will be able to send requests, via a Web interface of OPTEMIS system, by selecting the suitable satellites according to the type of the disaster in terms of sensor type and resolution (e.g. optical, radar, low, medium or high resolution). Then, our goal will be to make the most of these sensors to meet the needs of decision makers and search and rescue teams in the first days following the event. Eventually, an optimized acquisition tasking order will be particularly distributed to every collaborating satellite operators and/or space agencies. This work can be addressed as a multi-acquisitions by non-homogeneous satellites mission planning problem, where the inputs are requests defined by a polygonal region of interest (ROI), a set of satellites which will collaborate, a time horizon and a weight. The weights allow us, for instance, to give priority to urban areas over uninhabited ones. The model seeks to maximize the weighted-covered surface over the time horizons requested by the users while insuring the feasibility of the mission and a sufficient image quality. The non-homogeneous satellites and their different constraints (sensor type, field of view, satellite agility and maneuverability, and other particular satellite constraints) result in a complex combinatorial problem.

In this work, we suppose that each satellite operator considers the disaster requests as highest priority, and that the ROIs are relatively small. Therefore, it can be assumed that the constraints of on-board memory, energy, and ground station availability are satisfied. Other considerations concerning the fairness of resource sharing and the task distribution equity may be more appropriate for post-disaster monitoring of wide areas (non-urgent tasks).

For each satellite, we will need to predict its orbit passes over the requested areas and compute its observation opportunities based on its viewing characteristics. Then, we will try to plan acquisitions for each satellite by optimizing the whole mission with a metaheuristic optimization method. This process is suitable for the problem of co-constellation mission planning for disaster response, according to the general assumption for urgent disaster acquisitions that these tasks shall be put at the highest priority and shall override the existing plan in the acquisition queue.

II. Literature review

A lot of research has already been done on mission planning and scheduling for a single (agile) satellite and some other research papers deal with a cooperative satellite constellation. [START_REF] Sun | Multi-objective Mission Planning Problem of Agile Earth Observing Satellites[END_REF] proposed a multi-objective evolutionary algorithm to address four criteria: total number of acquisitions, total profit, image quality (minimizing the observation angles) and load balancing of satellites. In [START_REF] Bianchessi | Planning and scheduling algorithms for the COSMO-SkyMed constellation[END_REF], a constructive planning and scheduling algorithm with look-ahead and back-tracking abilities was developed for the COSMO-SkyMed constellation. The goal was to maximize the number of images taken and transmitted regarding downlink operational constraints. [START_REF] Bonnet | Multi-satellite Mission Planning Using a Self-Adaptive Multi-agent System[END_REF] split the request into a set of elementary areas that a satellite can acquire in a single shoot and solved the problem with a self-adaptive multi-agent system. In [START_REF] Li | A Decomposition-Based Algorithm for Imaging Satellites Scheduling Problem[END_REF], the problem was decomposed into a task assignment main problem and a single-satellite scheduling sub-problem. A Simulated Annealing heuristic was used for the single-satellite scheduling problem and the results were passed on an ant colony optimization algorithm to guide the assignment search process. A more comprehensive and comparative literature review can be found in [START_REF] Mitrovic-Minic | Collection Planning Management: Multi-Satellite Collection Scheduling[END_REF]. However, some of these works rely on a predefined set of strips without a clear indication on how they were computed, and so address only one part of the problem. In some papers, the authors split the area to cover in regular and contiguous strips just once, without considering the dynamic nature of the overlapping optimization problem over the days, which leaves no place for further optimization and can't be applied for a constellation of non-homogeneous satellites.

To the best of our knowledge, relatively few studies have investigated a computational model for a problem similar to ours. In [START_REF] Liu | Optimizing large area coverage from multiple satellite-sensors[END_REF], the authors used the circle intersection point set (CIPS) to derive a set of key representative geographic locations, similar to the facility location problem, based on the characteristics of the most restrictive sensor. This approach hugely reduces the combinatorial complexity; however, they rely on a preexisting module which gives the best acquisition opportunities for each interpolated point. Their optimization model also addressed weighted spectral and spatial resolution and off-nadir angles. They solve the problem as a set covering problem (SCP) with a branch-and-bound algorithm. In [START_REF] Liu | Satellite image collection modeling for large area hazard emergency response[END_REF], they extended the model to handle multiple polygonal requests with different weights. [START_REF] Niu | Multi-satellite scheduling approach for dynamic areal tasks triggered by emergent disasters[END_REF] proposed a dynamic segmenting algorithm for the multi-request problem for semi-agile satellites: for each orbit pass which can intersect with a requested area, they built a list of possible acquisitions by discretizing the roll angle depending on the sensor's characteristics. Then, they developed a greedy heuristic which also handles the dynamic arrival of new emergency tasks for a robust rescheduling. [START_REF] Wang | A model, a heuristic and a decision support system to solve the scheduling problem of an Earth observing satellite constellation[END_REF] proposed a model for the multi-request mission planning problem including download windows, which was evaluated for the environmental and disaster monitoring in China as part of the CRESDA program. Their model can be used for both spot and polynomial weighted requests. For agile satellites, they discretized the pitch angle and retain three time windows (-θ max , 0, +θ max ) whereas how they handled the continuous roll angle to generate the possible acquisitions for each orbit pass is not clearly explained. Also, their model didn't allow more than one strip's acquisition for a given request in each orbit pass, which can lead to a sub-optimal use of agility but reduces the complexity of the problem. The authors then developed a heuristic with good computational performances whose main idea "is to assign the most conflict-free observation time window to the most desired request, supplemented by limited backtracking". In [START_REF] Hodgson | Modeling the Potential Swath Coverage of Nadir and Off-Nadir Pointable Remote Sensing Satellite-Sensor Systems[END_REF], the authors addressed the geometrical problem of modeling the ground observation opportunities and demonstrated how a "temporal series of re-centered map projections" can be used to compute the strips. Other interesting papers are [START_REF] Li | Two-Archive2 Algorithm for Large-Scale Polygon Targets Observation Scheduling Problem[END_REF] and [START_REF] Rabideau | Scheduling Results for the THEMIS Observation Scheduling Tool[END_REF] which addressed the target discretization with a grid-based approach which reduces the calculation of overlapping areas by computing intersections of grid point sets. Finally, [START_REF] Liu | An adaptive large neighborhood search metaheuristic for agile satellite scheduling with time-dependent transition time[END_REF] gave some interesting insights on the neighborhood structure for agile satellite scheduling and defined six removal operators and three insertion operators.

III. Problem modeling

The mission planning problem we aim to address is complex due to the non-homogeneous orbits, the mix between non-agile and (semi-)agile satellites and their different kinematics. Therefore, we develop a simplified model which responds to the basic requirements of our problem.

A. Assumptions

To reduce the complexity of the problem, the following assumptions are made:

Requests and strips

For Earth observation missions, there are two main types of requests: spot or polygon. Usually, a spot can be acquired in a single shot whereas a polygon needs to be split into smaller observations (strips) which will be merged. In this paper, we consider polygon-shaped requests and we also handle non-convex and overlapping requests.

Strips are chosen parallel to the satellite's ground track to avoid costly maneuvering transition between successive shots which would result in fewer acquisitions ( [START_REF] Lemaître | Selecting and scheduling observations of agile satellites[END_REF]). As we support overlapping requests, strips are not attached to one specific request and can cover several requests.

Image acquisitions

We only consider full-strip acquisitions (the satellite will continuously acquire the strip from one edge of the polygon to the other one). It could mean that the pitch agility will be underexploited in some scenarios. One reason to avoid partial acquisitions of strips and the overuse of agility is that having too many different observations could result in a non-homogeneous quality and a lot of small, scattered acquisitions which are less valuable than continuous ones.

B. Input data: requests

We note the set of available satellites K = {k 1 , ..., k N K }. The input data is a set of requests R = {r 1 , ..., r N R }, each request r i consisting of:

• a polygonal region of interest (ROI) with lat/long coordinates: ROI Lat Lng (r i );

• a time interval for the mission: [T st art (r i ); T end (r i )];

• a set of sensors which will collaborate to perform the mission together: K(r i ) ⊆ K;

• a weight: w(r i ).

The problem will ultimately rely on geometric computations such as the intersection and the union of polygons, which are computationally expensive tasks. For efficiency, we decide to discretize each ROI with a grid-based approach ( [START_REF] Li | Two-Archive2 Algorithm for Large-Scale Polygon Targets Observation Scheduling Problem[END_REF]). We first project the polygonal ROI in the Lambert azimuthal equal-area space centered at the centroid of the ROI. Then, we compute a discrete set of points spaced by gridSpacing km covering the ROI. This parameter must be chosen carefully to reach a good balance between accuracy and performance, considering the size of the requests and the average swath width of the sensors (e.g. 10 km for very large requests with an average swath width of 25 km should be enough to compensate the side effects of discretization, whereas we would require 3 km for smallest ROIs).

Fig. 2 Discretization of the ROI

We note Points(r i ) the set of grid points for the request r i . N.B.: Even if two requests overlap, each request will still have its own set of points.

C. Strips modeling and computation

From these requests, we want to generate all the acquisition opportunities for the constellation. We note S(k) = {s k1 , ..., s k N S(k) } the set of strips for the satellite k. Each strip s ki will be modeled by: • A polygon with lat/long coordinates corresponding to the ground projection of the Earth surface it covers: Strip Lat Lng (s ki ); • A roll angle: φ(s ki ); • A duration: d(s ki );

• The discretized coverage of the ROIs, i.e. for each request r i , the set of the ROI's points it covers: ). The angles are discretized with a tunable step for each satellite (typically, 1°). Non-agile or semi-agile satellites will have only one time window (0°pitch angle).

∀r i ∈ R,
We define an acquisition as a strip s ki and one of its time window s ki [t]. For simplicity purposes, s ki [t] will be used to denote the acquisition, assuming it is also possible to retrieve the attached strip from this notation.

Orbit propagation and strip building

We describe below the algorithm used to build the strips for a given sensor. It is based on a fast propagation mode used to determine when a satellite intersects an ROI and a slow propagation mode used to build precise strips when an intersection is detected.

We first determine the set of requests the sensor is attached to and we merge their date intervals so that we know the useful timespan for propagation. We initialize the time increment of the propagation to 30 seconds.

1) We compute the lat/long ground position of the satellite by using an SPG4 library and its TLE. If the sensor is optical and is not in sunlight, go to the next date. 2) For each request attached to this sensor whose time horizon contains the current date:

• We compute the FOR of the sensor between the previous and the current position (at nadir, based on its maximum roll angle and its FOV angle) in the Lambert azimuthal equal-area euclidean space centered at the ROI's centroid. • We detect intersection between the FOR and the ROI in this space.

3) If the FOR intersects at least one request: if increment was 30 seconds, we decrease it to 0.5 seconds (or less if the ROIs are small) and go back to 1) at the previous date, else, for each discretized roll angle:

• We compute the "corridor" coordinates (boundaries of the strip at this position) based on the roll and the FOV angles in the Azimuthal equidistant space centered at this ground position ( [START_REF] Hodgson | Modeling the Potential Swath Coverage of Nadir and Off-Nadir Pointable Remote Sensing Satellite-Sensor Systems[END_REF]). • We project these coordinates back in lat/long. • For each intersecting request detected in 2):

-We project the lat/long corridor coordinates in the Lambert azimuthal equal-area space centered at the ROI's centroid. -We detect intersection in this space, if there is one: * If the strip was updated more than minGap seconds ago, we decide to build a new strip. First, we finalize the last strip (compute the discretized coverage of the ROIs it intersects and the list of its time windows) and then we create a new strip. * We extend the strip by adding the two corridor coordinates. 4) If the FOR intersects no request: if we were in slow mode, set increment back to 30 seconds. 5) Increment the date by increment and go back to 1). Figures 3, 4 and5 describe the propagation and strip building process. We use the minGap parameter to detect the end of a strip so as to decide whether we extend the strip or create a new one in case there was no intersection for some time. Indeed, for a given roll angle, if the last intersection comes from the previous orbit pass, from another non-overlapping request or from a previous position in a large non-convex ROI, we may decide it is better to switch off the sensor between two acquisitions.

Time windows

First, we compute the duration of the strip's acquisition based on the characteristics of the sensor. Then, we can compute the time windows by discretizing the pitch angle. For each possible pitch angle, we can approximate the starting date of the acquisition with basic geometry and compute the final pitch angle knowing the kinematics of the sensor.

Finally, we define the quality of an acquisition as a decreasing function of the roll and the initial pitch angles:

qualit y(s ki [t]) = 1 - |φ(s ki )| + |θ(s ki [t])| (|φ| + |θ|) max,K (1)
where (|φ| + |θ|) max,K is the worst distortion angle that could be produced by a sensor of the constellation. Thus, the quality lies between 0 (worst) and 1 (best). N.B.: We assume that requests are covered by sensors with similar resolution so that the value of the acquisition can be evaluated with its quality score. Otherwise, we should extend this quality function and consider the resolution.

Successive feasibility and strips clustering

For each satellite k, we assume that we have a function to get the minimum maneuvering duration between two successive acquisitions {strip s ki with time window s ki [t]} and {strip s k j with time window s k j [u]}. Then the physical feasibility constraint can be written:

f easibilit y(s ki [t], s k j [u]) := T st art (s ki [t]) + d(s ki ) + maneuvering k (s ki [t], s k j [u]) ≤ T st art (s k j [u]) (2)
Two strips are said to be in conflict if the feasibility constraint is violated for at least a pair of time windows (typically, the last time window of one strip and the earliest time window of the other strip). For each sensor, we cluster its strips so that we know for sure that there won't be any feasibility conflict between a strip from one cluster and a strip from another cluster. In most single-request scenarios with a straightforward geometry, a cluster would correspond to one orbit pass.

We note

C(k) = {C k1 , ..., C k N C(k)
} the set of clusters for the sensor k, with C kc ⊆ S(k).

D. Decision variables

For each satellite k, for each strip S(k), we can choose to plan zero (denoted by ∅) or one of its time windows, denoted by its index in TW(s). Therefore, we introduce the decision variables:

∀ k ∈ K, ∀ 1 ≤ i ≤ N S(k) , x ki = ∅ if strip s ki is not planned t if strip s ki is selected with time window s ki [t]

E. Objective

For each point p ∈ Points(r) belonging to the request r, we note Covering(p) the set of the planned acquisitionscoming from any sensor -that cover this point:

Covering(p) = k ∈K s k i ∈S(k), p ∈Points r (s k i ), x k i ∅ s ki [x ki ] (3)
We introduce the qualit yCoverage score of a point p as a function of the best quality that can be achieved by one of the planned acquisitions covering this point:

qualit yCoverage(p) =        0 if Covering(p) = ∅ f max acq ∈Covering(p) qualit y(acq) otherwise ( 4
)
where f is a tunable increasing function to balance the quality and the coverage, for instance a linear function with positive slope and y-intercept. Special cases are f = 1 where every covered point will have the same score regardless of the image quality and ∀q, f (q) = αq + β with β < 0 where worst quality leads to negative scores which will act as a soft constraint on the quality.

In this research, we are interested in short-term disaster response. The main consideration is maximizing the covered surface while optimizing the quality of the acquisitions, taking into account the weights of the requests. Therefore, we propose to maximize the following objective function:

max r ∈R w(r) × p ∈Points(r) qualit yCoverage(p) (5)
This objective function is designed to guarantee a first maximized coverage of the area in the first days following a major event. It is not intended to be used for monitoring as the qualit yCoverage score doesn't carry sufficient information for such purposes and cannot handle scoring for lower-quality and full-overlapping strips which won't improve the objective.

F. Constraints

First, we must ensure that the decision vector contains only valid time windows:

∀ k ∈ K, ∀ 1 ≤ i ≤ N S(k) , x ki ∈ {∅} ∪ 1, #TW(s ki ) (6)
The only constraint we consider is the intra-cluster feasibility in terms of maneuvering as described in 2:

∀ k ∈ K, ∀ C ∈ C(k), ∀(s ki , s k j ) ∈ C 2 , i j, x ki ∅, x k j ∅, T st art (s ki [x ki ]) ≤ T st art (s k j [x k j ]), ensure f easilibit y(s ki [x ki ], s k j [x k j ]) (7)
Other physical and operational constraints fall outside the scope of this work.

IV. Computational method

The mission planning problem involves a high combinatorial search space. To solve this problem, we propose below a heuristic to reach an approximate solution in polynomial time based on the Simulated Annealing algorithm.

A. Fast computation of the objective function

First, we explain how to make most of the model to compute quickly the objective function. For our problem, we can attach a list of quality scores to each point of each request and maintain the qualit yCoverage score of the point. Each time we decide to plan an acquisition, we loop over the points covered by the strip, we add the quality score of the acquisition to their list and update their qualit yCoverage score if needed. The reverse process happens when we decide to remove an acquisition. (Algorithm 1)

The objective function can now be computed efficiently for a state X by looping over the states of the points as described in Algorithm 2, with a constant-time Q C function.

B. Simulated Annealing

Simulated Annealing is a probabilistic metaheuristic motivated by the physical process of annealing in metallurgy. It involves heating a material to bring energy to the system and then cooling it down to a lower energy state. During the process, crystals are formed and their quality will depend on the thermodynamic properties of the material and the cooling rate. The Simulated Annealing heuristic (Algorithm 3) starts from a random state with high temperature and as the temperature is decreasing, it will probabilistically move to random neighboring states (little and localized changes in the decision vector) based on an acceptance criterion. The main idea is to allow deteriorating moves in the search space to escape from local optima, as long as the temperature (energy) is still high enough. While the energy is decreasing, such moves are less likely to be accepted and the heuristic will eventually favor improving states. The algorithm is stopped when the temperature has reached a lower bound or when the solution seems to be stabilized around a local optimum with no improvements over some iterations.

In this work, we use the Metropolis algorithm, inspired by the Boltzmann distribution, as the C : an improving move is always accepted and a deteriorating move is accepted with probability exp(-∆obj/T) where ∆obj is the magnitude of the deterioration between the current state and the proposed state and T is the current temperature. The temperature control will determine how much of the search space is explored and the quality of the solution. For the cooling scheme, we choose a geometric decrement α < 1 such that T ← αT at each iteration. In order to choose the initial temperature, we compute an estimate for ∆obj over successive moves by starting from random states. Then we reverse the acceptance criterion to get an initial probability of 0.70 to accept deteriorating moves.

Initialization

As our constraints are limited, we can start by selecting all the candidate strips of non-agile satellites with their single time window as there is no possible conflict between them. Then we randomly loop over the strips of agile sensors and try to plan each strip with a random time window while ensuring the feasibility (Algorithm 4). The C F (C, s ki [t]) function checks that the constraint is not violated if we plan s ki [t] belonging to cluster C.

Neighborhood structure

Metaheuristics such as Simulated Annealing rely on the generation of neighboring states with little and localized changes which are easily reversible. For instance, it can be inserting a new acquisition in the planning or a limited sequence of deletions of planned acquisitions followed by an insertion. We use different random operators and give a simple overview of such operations in Algorithm 5.

Finally, the B function is just responsible to cancel the changes made during the neighbor's generation.

The algorithms were implemented in Java and run on a 2 GHz CPU. From pre-processing to post-processing, the above instances were individually solved in about 30 seconds with strong convergence requirements on the metaheuristic. (a) f (q) = 1, 99.7% coverage (b) f (q) = q, 94.8% coverage (c) f (q) = 3q -2, 62.3% coverage

VI. Conclusion

In this article, we proposed a mission planning methodology that allows the user to maximize the use of satellite resources for disaster response. The proposed approach optimizes the mission planning for a constellation of nonhomogeneous satellites. First, a mathematical model was introduced, supporting multiple requests with different weights and both convex or non-convex ROIs. Then, a simplified strip-building algorithm was proposed. This pre-processing part must rely on a good knowledge of the physical characteristics of each collaborating satellite (e.g. orbital elements, camera characteristic and maneuvering capability for feasibility assessment). Finally, the mission planning problem was solved efficiently with the Simulated Annealing metaheuristic. Results from numerical experiments provide precious insight into possible cooperation in the very first days following a major event.

Future work should address common operational constraints (storage capacity, energy, downlink availability...) to ensure the feasibility of the mission and data delivering on time. The model could also be extended to other types of missions like post-disaster monitoring by seeking fairness in resource sharing.

Fig. 1

 1 Fig. 1 Illustration of the different types of agility. Non-agile satellites can have a fixed roll offset, meaning that their observations are not centered on their ground track: satellite 2 has a left-side viewing ability. The pitch agility of satellite 4 allows two observations in one orbit pass.

Fig. 3 4Fig. 5

 35 Fig. 3 Orbit propagation: acquisition opportunities are identified by intersecting the FOR and the requests between two successive positions. As soon as there is one intersection, we go back to the previous position and start again in slow mode.

Fig. 6

 6 Fig. 6 Strips for 3 satellites, over 3 days (Bangladesh)

Fig. 7 Fig. 8 Fig. 9

 789 Fig. 7 Selected acquisitions for 3 different f functions balancing coverage and quality (Bangladesh)

  Points r i (s ki ) ⊆ Points(r i ); • A list of time windows TW(s ki ), each time window s ki [t] is defined by: -A starting date: T st art (s

ki [t]); -A pitch angle: θ(s ki [t]); -A final pitch angle: θ f (s ki [t]); -A quality score: qualit y(s ki

[t]

VII. Acknowledgement

The authors would like to thank the OPTEMIS team for their contributions and GISTDA for sponsoring this project.

C. Post-processing

Once the optimization is done, the solution may include some acquisitions which fully overlap with acquisitions of better quality. These acquisitions have no effect on the objective function (as the maximum quality is already reached for the points they cover) and were chosen randomly, the reason why we decide to remove such acquisitions.

The solution may also be a local optimum which can still be improved. Especially in terms of time windows where some residual pitch angles could be reduced or even set to zero. Therefore, we decide to run a second optimization pass, initialized with the solution given by the first pass. We consider that the strips selection problem has been solved and our only concern is now to affect better time windows for the selected strips to improve the total average quality, with the same constraint of intra-cluster feasibility. We solve this problem with a similar heuristic and as the combinatorial search space is reduced, the complexity decreased and this post-processing process is not expensive.

V. Experimental results

The proposed model is evaluated on real requests issued by official agencies. The first scenario comes from the Department of Disaster Management of Bangladesh after flooding affected a rectangular 16000 km 2 area and 6 million people. The second request, defined by a 2800 km 2 circular area, was sent after flooding in several Vietnamese provinces. In both cases, quick disaster response is needed to find and give assistance to casualties and to assess damage to housing, infrastructures (roads, dams, bridges...) and agriculture.

For the simulation, three high-resolution optical sensors are collaborating over a maximum timespan of 3 days. THEOS and DUBAISAT-2 are agile satellites and maximum acceptable roll and pitch angles are set to ±30°in the experiment. VNREDSAT-1 is a semi-agile satellite with a field of regard of ±30°roll angle.

Figure 6 shows the strips computed over the 3-day timespan for Bangladesh, with a ±1°step for the discretized roll angles. In this simple scenario, the number of clusters happens to be the number of orbit passes over the ROI. At that point, we know for sure that we can plan at least 2 non-conflicting acquisitions for THEOS, 3 for DUBAISAT-2 and 1 for VNREDSAT-1. Pitch agility and time windows, which are not depicted, can allow several observations belonging to the same cluster and thus greatly increase the combinatorial complexity of the mission planning problem.

Figure 7 illustrates the results for three different f function used to balance the coverage and the quality and Figure 8 gives more information about the scheduled acquisitions. When quality is not taken into account and every observation is worth the same with f (q) = 1, the constellation can cover 99.7% of the ROI in 2 days thanks to the pitch agility of THEOS (7 a and 8 a). When f (q) = q, the scores are proportional to the quality and it can be seen than the pitch agility is no more beneficial and the coverage is reduced to 94.8% (7 b and 8 b). Finally, we set f (q) = 3q -2 so as to increase the weight of image quality and to introduce a soft constraint: indeed, only acquisitions with quality in [0.66, 1] will give a positive score and the bigger areas are now covered with a good quality (roll angles of 5°, -7°and 1°) to the cost of coverage, which drops to 62.3% (7 c and 8 c).

Figures 9 and10 show similar results for the Vietnamese request.