Scalable Real-Time Planning and Optimization Software Complex for the Purposes of Earth Remote Sensing

17th IFAC Workshop on Control Applications of Optimization 17th IFAC Workshop on Control Applications of Optimization Yekaterinburg, Russia, 2018 of online 17th IFAC Workshop onOctober Control 15-19, Applications Available at www.sciencedirect.com Yekaterinburg, Russia, 2018 of Optimization 17th IFAC Workshop onOctober Control 15-19, Applications Optimization Yekaterinburg, Russia, October 15-19, 2018 Yekaterinburg, Russia, October 15-19, 2018

ScienceDirect

IFAC PapersOnLine 51-32 (2018) 451–455

Scalable Real-Time Planning and Scalable Real-Time Planning and Scalable Real-Time Planning and Scalable Real-Time Planning and Optimization Software Complex for the Optimization Software Complex for the Optimization Software Complex for the Optimization Software Complex for the Purposes of Earth Remote Sensing Purposes of Earth Remote Sensing Purposes of Earth Remote Sensing Purposes of Earth Remote Sensing ∗ ∗ ∗

Vasiliy Vasiliy Sazonov Sazonov ∗∗∗ Sofia Sofia Sazonova Sazonova ∗∗ Ivan Ivan∗Samylovskiy Samylovskiy ∗∗ ∗ Vasiliy Sapelkin Sazonov ∗∗ Stanislav Sofia Sazonova Ivan∗Samylovskiy Arseniy Budzinskiy Oleg ∗ ∗ ∗ Arseniy Sapelkin Budzinskiy Oleg Bulatov Bulatov Vasiliy Sazonov Sofia Sazonova Ivan ∗ Stanislav ∗Samylovskiy ∗ ∗ ∗ ∗∗ Arseniy Sapelkin Stanislav Budzinskiy Oleg Bulatov Matvey Egorov Frolov ∗ Oleg Morozov ∗ Igor ∗∗ ∗ ∗ ∗ Matvey Egorov Morozov Frolov Arseniy Sapelkin Stanislav Budzinskiy Bulatov ∗ Oleg ∗ IgorOleg ∗∗ ∗∗ ∗∗ Matvey Egorov Morozov ∗ Tsvetkov Igor Frolov Vitaliy Lobanov ∗∗ Vjatcheslav ∗∗ ∗∗ ∗ Oleg Vitaliy Lobanov Vjatcheslav Tsvetkov Matvey Oleg Morozov ∗∗Tsvetkov Igor Frolov ∗∗ ∗∗ VitaliyEgorov Lobanov Vjatcheslav Alexey ∗∗ ∗∗ Alexey∗∗Strochenkin Strochenkin Vitaliy Lobanov Vjatcheslav Tsvetkov ∗∗ Alexey Strochenkin ∗∗ Alexey Strochenkin ∗ ∗ Lomonosov Moscow State University, Moscow, Russia, (e-mail: ∗ Lomonosov Moscow State University, Moscow, Russia, (e-mail: State University, Moscow, Russia, (e-mail: [email protected]) ∗ Lomonosov Moscow [email protected]) State University, Moscow, Russia, (e-mail: ∗∗ Lomonosov Moscow [email protected]) and ∗∗ S.P. Korolev Rocket and Space Space Corporation Corporation “Energia”, “Energia”, Korolev, Korolev, [email protected]) ∗∗ S.P. Korolev Rocket Corporation “Energia”, Korolev, Russia ∗∗ S.P. Korolev Rocket and Space Russia S.P. Korolev Rocket and Space Corporation “Energia”, Korolev, Russia Russia Abstract: Abstract: During During this this work, work, aa program program complex complex containing containing specialized specialized libraries, libraries, desktop desktop Abstract: During thisservices work, is a program complex containing specialized libraries, desktop applications and web developed for the purposes of Earth remote sensing and applications and web services is developedcomplex for the containing purposes ofspecialized Earth remote sensing and Abstract: During this work, a program libraries, desktop applicationsspace and shooting web services is lifecycle. developedComplex for the purposes of Earthand remote sensing anda managing orders allows operators remote clients managing space shooting orders lifecycle. Complex allows operators and remote clients applications and web services is developed for the purposes of Earth remote sensing andaa managing space shooting orders lifecycle. Complex allows operators and remote clients functionality to create order for getting Earth images in different modes w.r.t. target satellite functionality to create orderorders for getting Earth images allows in different modesand w.r.t. targetclients satellite managing space shooting lifecycle. Complex operators remote a functionality to create order technological for getting Earth imagesAs in adifferent modes w.r.t. target to satellite ballistic data shooting processes. part complex, libraries form ballistic data toand and shooting technological processes. As adifferent part of ofmodes complex, libraries to form functionality create order for getting Earth images in w.r.t. target satellite ballistic (in dataprescribed and shooting technological processes.including As a part of complex, libraries Software to form optimal sense) shooting flight tasks optimal (in sense)technological shooting programs, programs, including flight tasks generation. generation. Software ballistic dataprescribed and shooting processes. As a part of complex, libraries toterrain form optimal (in prescribed sense) shooting programs, including flight tasks generation. Software architecture allows users to add new tasks, new shooting modes, new satellites and architecture allows users to add new tasks, new shooting modes, new satellites andSoftware terrain optimal (in prescribed sense) shooting programs, including flight tasks generation. architecture allows users to add new tasks, new shooting modes, new satellites and terrain segment items for earth remote sensing orbital segment items for perpective perpective earthnew remote sensing orbital constellations. constellations. architecture allows users to add tasks, new shooting modes, new satellites and terrain segment items for perpective earth remote sensing orbital constellations. segment items(International for perpective earth remote sensing orbital constellations. © 2018, IFAC Federation of Automatic Control) Hosting by Elsevier Ltd. All rights reserved. Keywords: Keywords: Optimization Optimization Methods; Methods; Numerical Numerical Methods Methods for for Optimization; Optimization; Applications Applications in in Keywords: Optimization Methods; Numerical Science; MethodsEarth for Optimization; Applications in Economics, Management and Environmental Remote Sensing Economics, Management Methods; and Environmental Science; Earth Remote Sensing Keywords: Optimization Numerical Methods for Optimization; Applications in Economics, Management and Environmental Science; Earth Remote Sensing Economics, Management and Environmental Science; Earth Remote Sensing 1. INTRODUCTION INTRODUCTION In this this paper paper we we concentrate concentrate on on the the Planning Planning AWS AWS (Auto(Auto1. In 1. INTRODUCTION In this paper we concentrate on the Planning AWS (Automated Work Station) providing shooting routes planning mated Work Station) providing shooting routes planning 1. INTRODUCTION In this paper we concentrate on the Planning AWS (Automated Work Station) providing shooting routes planning functionality. Modern Earth Earth remote remote sensing sensing systems systems operates operates aa huge huge functionality. Modern mated Work Station) providing shooting routes planning Modern Earth remoteofsensing operates a huge functionality. amount of of snapshots differentsystems types. These These snapshots amount snapshots ofsensing different types. snapshots Modern Earth remote systems operates a huge functionality. amount of snapshots of different types. These snapshots are used to satisfy of scalable set of orders. Thus, specialare used of to snapshots satisfy of scalable set of orders. Thus, specialamount of different types. These snapshots are to satisfy of scalable set of orders. Thus, izedused software components are needed needed to cover cover totalspeciallifecy2. ROUTES ROUTES PLANNING PLANNING ized software components are to total lifecyare of used to satisfy of scalable set of orders. Thus, 2. ized software components are needed toand cover totalspeciallifecycle order from creation to satisfying corresponding 2. ROUTES PLANNING cle of order from creation to satisfying and corresponding izedofsoftware components are neededscalability toand cover total lifecy2. ROUTES PLANNING cle orderarchieving. from creation to satisfying corresponding snapshots In addition, addition, of modern modern Shooting routes routes are are presented presented as as aa number number of of flight flight tasks tasks snapshots archieving. In scalability of cle of order from creation to satisfying and corresponding Shooting snapshots archieving. addition, scalability of leads modern remote sensing sensing systemsInw.r.t. w.r.t. satellites number to Shooting routes are presented as astorage numberand of flight tasks arrays (FTA) sending to satellite containing remote systems satellites number leads to snapshots archieving. In addition, scalability of modern (FTA) sending to satellite and containing remote systems w.r.t. number leads to arrays Shooting routes are needed presented as astorage number of actions flight tasks necessitysensing of these these systems partsatellites or full full automation. automation. Durarrays (FTA) sending to satellite storage and containing parameters that are to determine determine what and necessity of systems part or Durremote sensing systems w.r.t. satellites number leads to parameters that are needed to what actions and necessity of work, these systems part or fullcomplex automation. Durarrays (FTA) sending to satellite storage and containing ing current a scalable software to support parameters that are needed to determine what actions and at what times should be initiated to realize needed routes. ing current work, asystems scalablepart software complex to support necessity of these or full automation. Durwhat times should be initiated to realize needed routes. ing current work, a scalable software complex support at parameters that are needed to determine what actions and remote sensing orders lifecycle was created created andtoprepared prepared what times should be initiated to realize needed routes. remote sensing orders lifecycle was and ingoperation. current work, a scalable software complex support at remote sensing orders lifecycle was created andto prepared at what times should be initiated to realize needed routes. to Software complex consists of the structural to operation. complexwas consists of the remote sensingSoftware orders lifecycle created andstructural prepared to operation. Software complex consists of the structural parts demonstrated in Fig. Fig. as they they are presented presented in webwebparts demonstrated in 11 as are in to operation. Software complex consists of the structural parts demonstrated in Fig. 1 as they are presented in webservice for complex user. service for complex user. parts demonstrated in Fig. 1 as they are presented in webservice for complex user. service for complex user.

Fig. Fig. 1. 1. Program Program Complex Complex Web Web Service Service Main Main Screen Screen Fig. Fig. 2. 2. Observation Observation Program Program Structure Structure Fig. 2. Observation Program Structure Fig. 1. Program Complex Web Service Main Screen Fig. 2. Observation Program Structure Fig. 1. Program Complex Web Service Main Screen 2405-8963 © 2018, IFAC (International Federation of Automatic Control) Hosting by Elsevier Ltd. All rights reserved.

Copyright © under 2018 IFAC 451 Control. Peer review responsibility of International Federation of Automatic Copyright © 2018 IFAC 451 Copyright © 2018 IFAC 451 10.1016/j.ifacol.2018.11.426 Copyright © 2018 IFAC 451

IFAC CAO 2018 452 Vasiliy Sazonov et al. / IFAC PapersOnLine 51-32 (2018) 451–455 Yekaterinburg, Russia, October 15-19, 2018

Flight task arrays generate observation program (see Fig. 2) which is updated w.r.t. adding new orders, changing meteo information, orbit corrections etc. Corresponding scheme is presented in Fig. 3.

Fig. 5. Visibility stripe decomposed to a number of shooting stripes (1) Split the visibility along the path into the survey strips; (2) Determining the intersection of order areas with shooting bands; (3) Sweeping out configurations that do not meet the conditions of the order (shooting angle, sun angle, etc.)

Fig. 3. Observation Program and Related Structures

Output data: The list of stripes (shooting configurations) that can potentially be shooted within the specified time interval. Each strip is characterized by: (1) (2) (3) (4)

Fig. 4. Observation Program Information in Web Interface

Shooting start and end time Area Yaw, pitch and roll Continuity (strip is continuous if it is possible to switch to a stripe with higher priority without completing current stripe shooting) It is necessary for shooting stereopairs, stereotypes. (5) List of orders with regions intersecting with this strip: (a) Order priority Pi (b) The ratio Ki of order-stripe intersection area to toltal strip area

2.1 Shooting Types • Single snapshot; • Stepeo pair (two snapshots with corresponding angles); • Stepeo triplet (three snapshots with corresponding angles); • Single strip of prescribed square (parallel to flight track); • Square shooting (four stripes with corresponding angles); • Corridor (geodesic-based strip characterized by azimuth and length). In addition, planning library provide an opportunity to plan curved-based corridors based on splines of corrsponding orders.

Fig. 6. Resulting Stripes 2.3 Routes Planning Algorithm. Step 2 – Create Graph Corresponding to Obtained Configurations Graph Nodes For each strip obtained on the previous step, we create graph node characterised by weight parameter as follows:

2.2 Routes Planning Algorithm. Step 1 – Prepare Input Data

W =S

N 

10Pi Ki ,

(1)

i=1

Input data:

where

• S is stripe square, • Pi is priority of i−th order intersecting with stripe, Si • Ki = , where Si is square of intersection between S stripe and order polygon.

• data array of orders; • satellite trajectory.

Algorithm:

452

IFAC CAO 2018 Vasiliy Sazonov et al. / IFAC PapersOnLine 51-32 (2018) 451–455 Yekaterinburg, Russia, October 15-19, 2018

Graph Edges (1) Creation of all possible edges. The edge goes from the vertex V1 to the vertex V2 , if the beginning of V2 is later than the time of the beginning of V1 . (2) Edges weights. The edge weight Wedge is calculated based on the conditions described below. If the weight of the edge is negative, then such an edge is not included in the graph. dt is the time difference between the beginning of V2 and the end of V1 , measured in milliseconds. dm - the number of milliseconds required to rebuild the device from the V1 configuration to the V2 configuration (a) If dt > 0, dm < dt, then the device has time to completely remove the strip V1 , then it will be rebuilt on the strip V2 and begin to shoot V2 . Thus, from rebuilding, we win the weight of vertex V2 entirely. Hence, Wedge = V2 .

453

of the vertex, which is considered by the formula (1). T is the duration of the claimed V-band t survey. Wedge = V2 − V T

Fig. 9. . (3) Adding an ”initial” and ”final vertex” (a) A vertex A is added, which is connected to all other vertices by edges with weights equal to the weights of these vertices. This means that we can start shooting from any band without loss. (b) A vertex B is added. Each vertex is connected to it by an edge with a weight of 0. This means that after shooting any strip, we can go to the ”rest state” and do nothing more. .

Fig. 7. dt > 0, dm < dt (b) If dt > 0, dm > dt, then the device will not be able to reconfigure to the V2 configuration to start shooting at the specified time. Accordingly, he needs an extra time t = dm − dt. You can get this extra time either by completing the v1 strip earlier, or by starting to remove the V2 strip later. We will “cut” the lane with less weight, if both of them can be “interrupted”. If one of the bands is “inseparable”, then the second one is automatically selected. If both are inseparable, then we assume that there is no such edge: V = min (V1 , V2 ) , the comparison is based on the node weight (1). T is the duration of the claimed Vband survey. t Wedge = V2 − V , T i.e. when we perform the rebuilding V1 − V2 , we ”acquire” the value V2 for the target functional, but ”lose” in the rebuilding.

Finding Optimal Route After performing the actions described above, we got an oriented graph without cycles. In such a graph one can find the path of maximum length. To do this you need: (1) Perform topological sorting. Method of topological sorting with the help of a bypass in depth. A - the first vertex after sorting, B - the last one; (2) For each vertex v in topological sorting, we calculate the length of the longest path terminating at the vertex v by viewing the incoming arcs from the neighbors and adding the weight of the connecting edge to the maximum length in the records of these neighbors. If v has no incoming arcs, we assign the length of the longest path ending in v to zero; (3) The longest path to the vertex B is the desired sequence of surveys 3. CORRIDOR SHOOTING 3.1 Algorithm to find corridor shooting route Input data: Polygon points array pi = (lati , loni , alti )

Fig. 8. dt > 0, dm > dt (c) If dt < 0, then the bands intersect in time. In this case, if dm milliseconds is required to rebuild, then we ”lose” the time t = dm + |dt|. (d) We will “cut” the lane with less weight, if both of them can be ”interrupted”. If one of the bands is “inseparable”, then the second one is automatically selected. If both are inseparable, then we assume that there is no such edge: V = min V1 , V2 , the comparison is based on the weight 453

Fig. 10. Initial Polygon Algorithm • Building linear skeleton of the polygon;

IFAC CAO 2018 454 Vasiliy Sazonov et al. / IFAC PapersOnLine 51-32 (2018) 451–455 Yekaterinburg, Russia, October 15-19, 2018

• Representation of the skeleton as the graph; • Finding graph path of maximal length.

in the right direction). For example, if we use “regular” corridores to cover curve track the result is as shown on Fig. 13.

Output data: Characteristic line points array pi = (lati , loni , alti ) 3.2 Building Sceleton As a first step, we build sceleton according to Felkel and Obdrzalek (1998) algotithm: (1) For each vertex, finding intersections of the bisector with the bisectors of neighbour vertices; (2) First point, where the edge is contracted, is added to skeleton; (3) Connect this points with vertices of contracted edge; (4) Rebuild polygon, moving vertices along the bisectors; (5) Loop: go to step (1) until possible. Output data (see Fig. 11):

Fig. 13. Covering curve line with a collection of “regular” corridores

• Skeleton vertices array; • Skeleton edges array.

Planning library provides also coverage shown on Fig. 14.

3.3 Finding Characteristics Line

Fig. 11. Constructed polygon sceleton • Characteristic line (see Fig. 12) is closest path between two furthest skeleton points; • Bellman-Ford algorithm Bellman (1958) is used for finding path between each two points; • Length is the measurement of edges.

Fig. 14. Covering curve line with a collection of “curvebased” corridores The structure of route parameters allows us to use 4 coefficients to define corridors: 2 for latitude and 2 for longitude. Hence, second order algebraic curves are used. 4.1 General Approach (1) (2) (3) (4)

Fig. 12. Constructed characteristic line 4. CURVE-STYLED CORRIDOR PLANNING

Construct polygon’s characteristic line; Subdivide it into smaller parts; Approximate each part with a second order curve; Construct corridors for each part.

4.2 Characteristic line subdivision

Modern earth remote sensing systems (see Ple (2012)) support cyrve-styled corridores to shoot rivers, state borders, roads etc. The goal is to cover the given thin polygon with a corridor that adapts to track the form of the polygon (i.e. curves 454

Input: Polygons characteristic line as a list of (latitude, longitude) pairs. We break the characteristic line into parts according to the curvature sign: each part will have either positive or negative curvature.

IFAC CAO 2018 Vasiliy Sazonov et al. / IFAC PapersOnLine 51-32 (2018) 451–455 Yekaterinburg, Russia, October 15-19, 2018

κ=

x y  − y  x

(x2

+

3/2 y 2 )

s(t) = s0 t + s1 t2 + s2 t3 .

(2)

Here, x and y are latitude and longitude as a functions of length from the rouute start point, respective; derivatives are taken with respect to arc length Output: Subdivision of the characteristic line as a list of lists of (latitude, longitude) pairs An example of subdivision is shown on Fig. 15, Fig. 16.

455

(4)

To fit a second order curve to a part of the polygons characteristic line with constant sign curvature, we find the point of maximum curvature. The curve is then uniquely defined by 3 points: first point of the curve, maximum curvature point, last point of the curve. 4.4 Corridor Construction Input: • Starting time; • List of coefficients, 4 for each characteristic line part

Output: List of polygons that represent corridors

The actual corridor is constructed so that the satellites point of view follows the second order curve. To this end, we need to find the relation between arc length of the curve and time.

Fig. 15. Original polygon for corridor shooting - curvature sign changed

Fig. 17. Corridors created to cover polygons with changing and constant curvature sign 5. CONCLUSION Fig. 16. Polygon with changing curvature sign separated into two corridores 4.3 Approximation with a second order curve Input: Part of the characteristic line with constant sign curvature as a list of (latitude, longitude) pairs Output: 4 coefficients that define a two-dimensional second order curve We are looking for an approximation of the form x(s) = x0 + x1 s + z2 s2 , y(s) = y0 + y1 s + y2 s2 ,

A collection of libraries, desktop applications and webservises is developed to form a complex that allow operator a functionality to create order for getting Earth images in different modes w.r.t. target satellite ballistic data and shooting technological processes. As a part of complex, libraries to form optimal (in prescribed sense) shooting programs, including flight tasks generation. Software architecture allows users to add new tasks, new shooting modes, new satellites and terrain segment items for perpective earth remote sensing orbital constellations. REFERENCES

(3)

where x(s), y(s) are latitude and longitude, s is route length, x0 and y0 are curve start point initial latitude and longitude, x1 , y1 , x2 , y2 are the spline coeffitients to be found. In addition, we use the following interpolation for route length w.r.t. time t : 455

(2012). Pleiades Imagery User Guide. Astrium GEOInformation Services. Bellman, R. (1958). On a routing problem. Quarterly of Applied Mathematics, 16(1), 87–90. Felkel, P. and Obdrzalek, S. (1998). Straight skeleton implementation. SCCG 98: Proceedings of the 14th Spring Conference on Computer Graphics, 210–218.