User Guide

Running temporal logic control synthesis is more than just a single line of code. The user needs to

  • specify the system dynamics, which falls into the following categories:

    • continuous/discrete-time dynamical systems without controls \(\dot{x}=f(x)\)/\(x^+=f(x)\);

    • continuous/discrete-time switched systems in the form of \(\dot{x}=f_u(x)\)/\(x^+=f_u(x)\), where the control input \(u\in\mathbb{Z}\) indicates the system mode;

    • continuous/discrete-time control systems in the form of \(\dot{x}=f(x,u)\)/\(x^+=f(x,u)\);

  • specify the control specification, which can be

    • invariance: \(\varphi=\Box a\), or

    • reachability: \(\varphi=\lozenge a\), or

    • an LTL formula that can be translated into a DBA,

  • save the resulting control synthesis strategies in either

    • the MATLAB file format or

    • the HDF5 file format

  • prepare a simulation file to examine the closed-loop performance.

This manual will walk you through a typical formal control synthesis procedure, which includes the above points.

Prepare the Main Program

The first step to solving a control synthesis problem is to create a C++ main program that manages the workflow from the inputs of system dynamics and specifications to the output of control strategies. Such a file typically consists of:

  • read the control specification,

  • define the system dynamics,

  • assign labels to the areas in the state space by the labeling function,

  • call a proper control synthesis algorithm, and

  • save control strategies.

In the following sections, we will show you how to write each part.

Structure Your Specification

To perform control synthesis w.r.t. a DBA specification, a structured file describing the DBA is required currently. An acceptable specification file should contain

  • metadata about the DBA such as the equivalent LTL formula, the number of states, accepting states, atomic propositions and so on.

  • all the transitions (or edges) between states in the DBA, which are given in the order of from state, proposition of the edge, to state.

The format of the file is consistent with the Hanoi Omega Automata (HOA) format, which is a text-based format for many types of finite automata over infinite words, i.e. \(\omega\) -automata. To write such a specification file, the user can use the online translator Spot, which converts LTL formulas to the corresponding \(\omega\) -automata.

Consider the generalized Buechi specification \(\varphi_{\rm gb}=\Box\lozenge g_1\wedge\Box\lozenge g_2\) whose DBA translation is

_images/gb2.png

The DBA of the generalized Buechi formula.

The following structured text file is a transciption of the above DBA, which can be parsed by ROCS.

 1Name=GFg1 & GFg2
 2AP=g1(1) g2(2) others(0)
 3acc=0
 4nNodes=3
 5nAP=2
 6init=0
 7###(q_from),(prop),(q_to)###
 8M=[
 90,3,0
100,0,1
110,1,1
120,2,2
131,3,0
141,0,1
151,1,1
161,2,2
172,1,0
182,3,0
192,0,2
202,2,2
21]

Lines 1-6 are the metadata. The keywords (case sensitive) Name, AP, acc, nNodes, nAP and init should be used to indicate the corresponding LTL formula, the atomic proposition, the accepting states, the number of DBA states, the number of atomic propositions, and the initial state, respectively. The order of presenting these information does not matter.

Note

To use ROCS, there can be multiple accepting states while the initial state is unique. Separate different accepting states by comma in the specification file, e.g., acc=0,3,5.

Starting from line 8, transitions of the DBA are given in the form of a matrix enclosed by M=[]. Each line represents a transition from the starting state (column 1) to the ending state (column 3), as commented in line 7. The propositional formula enabling each transition is given in column 2 by the decimal value of its binary code.

The binary encoding of a propositional formula projects the formula into a binary code, which contains nAP bits (or binary digits). Each bit corresponds to an atomic proposition. The value of the bit is 1 when the atomic proposition is true; otherwise it is 0. In the above exmaple, \(g_1\) is represented in the lower bit and \(g_2\) the higher bit. Thus, \(g_1\wedge g_2=\set{11}=\set{3}\), and \(g_1=\set{10,11}=\set{2,3}\) and so forth.

In the main program, read the sepcification file by using read_spec:

1    std::cout << "Loading the specification...\n";
2    std::vector<rocs::UintSmall> acc;
3    std::vector<std::vector<rocs::UintSmall>> arrayM;
4    rocs::UintSmall nAP = 0, nNodes = 0, q0 = 0;
5    if (!rocs::read_spec(specfile, nNodes, nAP, q0, arrayM, acc)) 
6	std::exit(1);

Note

For the invariance and reachability control specifications, specialized control algorithms can be used directly. Specifically, instead of using a specification file that represents its corresponding DBA, the user can alternatively (recommended)

  • call methods invariance_control and/or reachability_control of the class CSolver if the specification-guided engine is used. In this case, the goal and avoid areas in the state space of the system needs to be specified firstly by using the init method.

  • call methods invariance and/or reachability of the class DSolver if the abstraction-based engine is used. In this case, the user needs to prepare a set of goal states (of type std::vector<size_t>) before using those methods.

Structure Your Dynamics

The system dynamics can be written either

  • directly in the main file or

  • in a separate header file.

Take the mobile robot problem in the Case Studies as an example, the vehicle dynamics is definded as a struct in the following header file car.hpp:

 1#ifndef _robotcar_h
 2#define _robotcar_h
 3
 4
 5#include <vector>
 6#include <cmath>
 7
 8
 9/* user defined dynamics */
10const double h = 0.3;  // sampling time
11struct carde {
12
13    static const int n = 3;  // system dimension
14    static const int m = 2;
15
16    /**
17     * Discrete-time dynamics
18     * @param h sampling time
19     * @param x system state: [x,y,theta], n=3
20     * @param u control array (size of 2, velocity and steering angle)
21     * @param nu the number of different control values
22     */
23    template<typename S>
24    carde(S &dx, const S &x, rocs::Rn u) {
25	double alpha, at, r;
26	alpha = atan(tan(u[1]) / 2);
27
28	if (std::fabs(u[0]) < 1e-6) {
29	    dx = x;
30	} else if (std::fabs(u[1]) < 1e-6) {
31	    dx[0] = x[0] + u[0]* cos(x[2])*h;
32	    dx[1] = x[1] + u[0]* sin(x[2])*h;
33	    dx[2] = x[2];
34	} else {
35	    at = alpha + u[0] * tan(u[1]) * h / 2;
36	    r = 2 * sin(u[0]*tan(u[1])*h/2) / (cos(alpha) * tan(u[1]));
37
38	    dx[0] = x[0] + r * cos(at + x[2]);
39	    dx[1] = x[1] + r * sin(at + x[2]);
40	    dx[2] = x[2] + u[0] * tan(u[1]) * h;
41	}
42    }
43
44}; // struct carde
45
46#endif

Note

Define your dynamics by using the typename S as in the example.

In the main file, create a system object by using the carode struct:

setting up a discrete-time system
 1    const double theta = 3.5;  // M_PI
 2    double xlb[] = {0, 0, -theta};
 3    double xub[] = {10, 10, theta};
 4    /* Control values */
 5    double ulb[] = {-1.0, -1.0};
 6    double uub[] = {1.0, 1.0};
 7    double mu[] = {0.3, 0.3};
 8    /* Control system */
 9    rocs::DTCntlSys<carde> carvf("DBA", h, carde::n, carde::m);
10    carvf.init_workspace(xlb, xub);
11    carvf.init_inputset(mu, ulb, uub);

The mobile robot is a discrete-time control system, and therefore the system type DTCntlSys is used in line 9. Please use the correct type for your system:

  • Discrete-time control systems: DTCntlSys;

  • Discrete-time switched systems: DTSwSys;

  • Discrete-time systems without controls: DTSys;

  • Continuous-time control systems: CTCntlSys;

  • Continuous-time switched systems: CTSwSys;

  • Continuous-time systems without controls: CTSys.

Typically, line 10 and 11, which initialize the state space and the set of sampled control inputs, respectively, are followed immediately after the control system object is created. For systems without controls or switched systems, only line 10 needs to be applied.

To handle continuous-time systems, high-order Taylor expansion is used to compute a validated reachable set after one sampling time for each cell in the system state partition, which is required by the control synthesis algorithms. The class flowTaylor is designed to compute validated reachable sets. Here is an example code for setting up a continuous-time system.

setting up a continuous-time system
 1    /* set the sampling time and disturbance */
 2    double tau = 0.05;
 3    double delta = 0.01;
 4    /* parameters for computing the flow */
 5    int kmax = 5;
 6    double tol = 0.01;
 7    double alpha = 0.5;
 8    double beta = 2;
 9    rocs::params controlparams(kmax, tol, alpha, beta);
10    
11    /* define the control system */
12    rocs::CTSys<vdpode> vdproa("VanDerPol", tau, vdpode::n, delta, &controlparams);
13    vdproa.init_workspace(xlb, xub);
14    vdproa.allocate_flows();

Compared to setting up a discrete-time system, there are extra lines:

  • Line 14 allocates memory for an flowTaylor object for validated reachable set computation. Since new memory is allocated, the following line of code shall be used at the end of the main program.

vdproa.release_flows();

  • Line 5-9 provide the parameters for validated reachable set computation (see C++ APIs and [LL2020] for the explanation of the parameters.)

Note

The parameters used in setting up a continuous-time system work for many dynamical systems. You may need minor modifications, but it is not suggested to change them without using them for your system first.

Define the Labeling Function

A labeling function connects the atomic propositions or propositional formulas given in the control specification to certain areas in the state space of the dynamical system. One way to define a labeling function is through the use of a lambda function.

If the specification-guided engine is to be used, then, first of all, we need to create nNodes CSolver objects, where nNodes is the number of states in the DBA specification. This is due to the structure of the overall DBA control strategy: the DBA is embedded to update the controller memory of the current DBA state. The implementation of such a control strategy in the specification-guided engine works on nNodes copies of CSolver objects. Each such object holds and controls the partition of a copy of the state space that corresponds to one of the DBA states. The control strategy for each DBA state will then be computed and stored in the corresponding CSolver object, which can be applied whenever its corresponding DBA state is activated.

As a result, the labeling procedure needs to performed for all the CSolver objects by the member function labeling in the solver class CSolver. Here is an example with the same mobile robot model:

Labeling with CSolver
 1    /* Obstacles */
 2    rocs::UintSmall nA = 4;
 3    double olb[4][3] = {{1.6, 4.0, -theta},
 4			{3.0, 5.0, -theta},
 5			{4.3, 1.8, -theta},
 6			{5.7, 1.8, -theta}};
 7    double oub[4][3] = {{5.7, 5.0, theta},
 8			{5.0, 8.0, theta},
 9			{5.7, 4.0, theta},
10			{8.5, 2.5, theta}};
11    /* Goals */
12    rocs::UintSmall nG = 3; // # of goals
13    double glb[3][3]= {{1.0, 0.5, -theta},
14		       {0.5, 7.5, -theta},
15		       {7.1, 4.6, -theta}};
16    double gub[3][3]= {{2.0, 2.0, theta},
17		       {2.5, 8.5, theta},
18		       {9.1, 6.4, theta}};
19    /* Initialize the set of S-domains */
20    std::vector<rocs::CSolver*> w(nNodes);
21    std::vector<rocs::SPtree*> sdoms(nNodes);
22    rocs::UintSmall labels[]{4, 2, 1}; // corresponding to goal1,2,3.
23    auto init_w = [&carvf, &nProps, &labels, &arrayM,
24		   &olb, &oub, &nA,
25		   &glb, &gub, &nG](std::vector<rocs::CSolver*> &w,
26				    rocs::UintSmall i,
27				    rocs::UintSmall oid[]) {
28		      w[i] = new rocs::CSolver(&carvf, nProps, rocs::RELMAX);
29		      for (rocs::UintSmall j = 0; j < nA; ++j)
30			  w[i]->init(rocs::AVOID, olb[j], oub[j]); // avoid areas
31		      for (rocs::UintSmall j = 0; j < nG; ++j)
32			  w[i]->labeling(glb[j], gub[j], labels[j]); // labeled areas
33		      w[i]->set_M(arrayM[oid[i]]);
34		  };

Line 23-34 is a lambda function that creates nNodes copies of CSolver objects as well as assigns the labels to specify goal and avoid areas in each of the object. Specifically,

  • line 1-18 specify the 3d rectangular goal and avoid areas.

  • line 22 specifies the labels of all the areas involved,

  • Line 28 creates nNodes CSolver objects, which are labeled in line 30, 32, and

  • the pointers to all these CSolver objects are collected in a vector w in line 20.

Similarly, if the abstraction-based engine is to be used, then an abstraction of the system needs to be defined first, and a lambda function can be used for labeling. In the following code example,

  • the first 2 lines define and initialize an abstraction of the mobile robot system,

  • line 9 constructs all the transitions between abstraction states,

  • line 18-35 assign different labels 4, 2, 1 to 3 rectangular goal areas, and the lambda function defined in line 22-34 serves as the labeling function.

 1    rocs::abstraction<rocs::DTCntlSys<carde>> abst(&car);
 2    abst.init_state(eta, xlb, xub);
 3    std::cout << "The number of abstraction states: " << abst._x._nv << '\n';
 4    /* Compute abstraction */
 5    /* Robustness margins */
 6    double e1[] = {0,0,0};
 7    double e2[] = {0,0,0};
 8    tb = clock();
 9    abst.assign_transitions(e1, e2);
10    // abst.assign_transitions();
11    te = clock();
12    float tabst = (float)(te - tb)/CLOCKS_PER_SEC;
13    std::cout << "Time of computing abstraction: " << tabst << '\n';
14    std::cout << "# of transitions: " << abst._ts._ntrans << '\n';
15    /*
16     * Assign labels to states: has to be consistent with the dba file.
17     */
18    double goal[3][4] = {{1.0, 2.0, 0.5, 2.0},
19			 {0.5, 2.5, 7.5, 8.5},
20			 {7.1, 9.1, 4.6, 6.4}};
21    rocs::UintSmall nGoal = 3;
22    auto label_target = [&goal, &nGoal, &abst, &eta](size_t i) {
23		      std::vector<double> x(abst._x._dim);
24		      abst._x.id_to_val(x, i);
25		      double c1= eta[0]/2.0; //+1e-10;
26		      double c2= eta[1]/2.0; //+1e-10;
27		      boost::dynamic_bitset<> label(nGoal, false); // n is the number of goals
28		      for(rocs::UintSmall i = 0; i < nGoal; ++i) {
29			  label[2-i] = (goal[i][0] <= (x[0]-c1) && (x[0]+c1) <= goal[i][1] &&
30				      goal[i][2] <= (x[1]-c2) && (x[1]+c2) <= goal[i][3])
31			      ? true: false;
32		      }
33		      return label.to_ulong();
34		  };
35    abst.assign_labels(label_target);
36    std::cout << "Specification assignment is done.\n";

Note

The labels used in the labeling function should be obtained by the same binary encoding as the ones in the specification.

Perform Synthesis

To perform the control synthesis algorihm with the specification-guided engine, the method dba_control should be used. Let us still look at the mobile robot example. The following line of code will perform the control synthesis w.r.t. the user-defined DBA specification for the mobile robot system.

1    rocs::dba_control< rocs::DTCntlSys<carde> >(w, &carvf, sdoms, nNodes,
2						isacc, init_w, oid, e);

The first input argument w is a vector of pointers to the CSolver objects, which has been created in Labeling with CSolver when we define the labeling function. Each such object is associated to a different DBA state. In some cases where only a sub-graph of the entire DBA is concerned with, only the pointers to the CSolver objects related to the sub-graph need to be collected in w. Hence, the size of w could be smaller than nNodes. After the control sythesis finishes, a control strategy for each DBA state will be stored in the corresponding CSolver object. In this way, w is also a returning argument from which we can save the control strategies for the corresponding CSolver objects to a file.

The input argument sdoms is a vector of pointers to the system state spaces possessed by all the CSolver objects, and the size of sdoms is nNodes. Since non-uniform partitions will be generated during control synthesis, a tree structure is used for representing each of the state space. Winning set will be marked in such a tree data structure. The sdoms is needed even when only a sub-graph of the DBA is considered, and an index mapping oid will be used to convert the index of the vector w to the corresponding index in sdoms.

The input function init_w is the lambda function we defined earlier for labeling, and e is the partition precision specified by the user. The argument isacc is a vector of binaries of size nNodes, which marks the accepting states. It can be obtained by using the acc after reading the specification file. Here is a code example:

1    boost::dynamic_bitset<> isacc(nNodes, false);
2    for (rocs::UintSmall i = 0; i < acc.size(); ++i)
3	isacc[acc[i]] = true;

To use the abstraction-based engine, there will be 4 steps:

  • construct a DBA from the information read from the specification file,

  • load the abstraction that has been constructed,

  • generate the product system of the DBA and the abstraction, and

  • run a Buechi game algorithm over the product.

Here is an example code (for the mobile robot model):

 1    /**
 2     * Solve a Buchi game on the product of NTS and DBA.
 3     */
 4    std::cout << "Start solving a Buchi game on the product of the abstraction and DBA...\n";
 5    rocs::BSolver solver;
 6    solver.construct_dba((int)nAP, (int)nNodes, (int)q0, acc, arrayM);
 7    tb = clock();
 8    solver.load_abstraction(abst);
 9    solver.generate_product(abst);
10    solver.solve_buchigame_on_product();
11    te = clock();
12    float tsyn = (float)(te - tb)/CLOCKS_PER_SEC;
13    std::cout << "Time of synthesizing controller: " << tsyn << '\n';

Save control synthesis results

In order to run closed-loop simulations, you need to write several lines of code to save control synthesis results as well as problem settings to one of the following two data formats:

  • the HDF5 format (a .h5 file), and

  • the MATLAB data format (a .mat file).

It is encouraged to use the HDF5 format although .h5 and .mat are interchangeable. This is because HDF5 is open distributed and also the base for the .mat format. For non-MATLAB users, the installation of HDF5 is necessary in order to use the Python interfaces. Besides, the generation of .mat files requires the installation of MATLAB.

To save the control strategies that are obtained by using the specification-guided engine into files,

  • use the method write_csolvers_to_h5 for writing .h5 files, or

  • use the method write_results_to_mat for writing .mat files.

The control strategy for each DBA state (or associated with each CSolver object) will be written into a single file. In other words, if there are \(|Q|\) states in the DBA, then \(|Q|\) number of .h5 files will be generated.

To save problem settings (such as system states, discretization precision) as well as the control strategies that are obtained by using the abstraction-based engine into .h5 files, the user can refer to the following example code:

1    rocs::h5FileHandler ctlrWtr(datafile, H5F_ACC_TRUNC);
2    ctlrWtr.write_problem_setting< rocs::DTCntlSys<carde> >(car);
3    ctlrWtr.write_array<double>(eta, carde::n, "eta");
4    ctlrWtr.write_2d_array<double>(abst._x._data, "xgrid");
5    ctlrWtr.write_discrete_controller(&(solver._sol));

Note

Writing to .mat files is not supported for the control systhesis with the abstraction-based engine.

Compile and run

So far, we have had the main program ready for control synthesis. The last step before running the actual control synthesis algorithms is to compile the main program. For this purpose, we suggest to prepare a Makefile under the same folder as the main program. The user can copy the Makefile from any existing examples and modify the compiling rules as needed. For example, if you have a main program called main.cpp, then modify the rule to

main: main.o $(OBJS)
        $(CXX) $(CXXFLAGS) $(INCS) $^ -o $@ $(LDFLAGS) $(LIBS)

Then type make main in the terminal.

Simulate Your Result

It is important to run closed-loop simulation to test the control performace once you have the control synthesis results generated. In this last section, we will see how to interpret the data saved in .h5 or .mat files and how to use them to do simulation.

You may expect the following data variables in each .h5 or .mat file generated by using the specification-guided engine:

  • U: a 2d array representing the set of sampled control inputs. Each row is a control input value, and the number of rows is the number of control inputs (let’s say \(|U|\)). The input dimension determines the number of the columns.

  • X: an \(n\times 2\) matrix representing the state space, where \(n\) is the system dimension.

  • ts: the sampling time.

  • pavings: an array of intervals representing the non-uniform partition of the system state space. Each row represents an interval (or cell). There are \(2n\) columns, which specify the lower and upper bounds in each dimension. The number of rows is the total number of cells (let’s say \(N\)) in the partition.

  • tag: a binary array of size \(N\) indicating the winning set of the control problem. If a cell belongs to the winning set, the corresponding row is 1; otherwise 0.

  • ctlr: a 2d binary feedback control table with \(N\) rows and \(|U|\) columns. The value 1 at row \(i\) column \(j\) means that the \(j\) th control input can be used to realize the given control specification when the current system state belongs to the \(i\) th interval (or cell).

Here we show how to use the above data for simulation by a python code example:

 1if(any(ctlr[q][x_id, :])):
 2     uset = np.argwhere(ctlr[q][x_id, :]).squeeze()  # get the indices of valid input
 3else:
 4     print("No valid control input.")
 5
 6if(uset.size > 1):
 7     uid = rng.choice(uset, 1)  # randomly pick one
 8else:
 9     uid = int(uset)
10u = U[uid, :].squeeze()

Since there is one .h5 or .mat file for each DBA state, all the control tables ctlr are put into a Python list after reading from the files. In line 1 of the above code, ctlr[q] is the control table corresponds to the DBA state \(q\).

In the .h5 file generated by using the abstraction-base engine, the data variables U, X, and ts are also included. Besides, there are the following additional data variables:

  • eta: the discretization precision with the same dimension as the system state.

  • xgrid: an array of discretized states of the system. Each row represents a state. There are \(N\) rows and \(n\) columns, where \(N\) is the number of the discrete states and \(n\) is the system state dimension.

  • WinSet: the winning state indices. The size of WinSet is the size \(N_{\rm win}\) of the winning set.

  • OptCtlr: an array of control pairs \((q,u)\), where \(q\) and \(u\) are the indices of DBA state and control input, respectively. To look for a control input \(u\) to realize the specification, we need the current system state and DBA state.

  • nts_ctrlr: a table with \(N_{\rm win}\) rows and 3 columns. The \(i\) th row stores the keys for looking up the corresponding control input for the \(i\) th state in the winning set WinSet. The columns represent the number of rows in OptCtlr to be searched, the label of the current system state, and the starting row in OptCtlr, respectively.

  • encode3: an index mapping that maps the indices of the states in xgirds to the indices in nts_ctrlr, i.e., matches row indices between xgrids abd nts_ctrlr.

  • q_prime: an array that encodes the DBA state update mechanism.

Let us assume that the current system state \(x\) corresponds to the \(i\) the state in xgrid. To determine the control input for xgrid[i], we should look at the encode3[i] th row of nts_ctrlr, which is p5=nts_ctrlr[encode[i],...]. Then p5[2] gives the starting row in OptCtlr to look for the actual control input, and p5[0] is the number of row we should search from the starting row. If the current DBA state index is \(q_c\), then we search among these rows to find the pairs with \(q=q_c\). The control input \(u\) paired with \(q_c\) is the control input we are looking for \(x\). Here is a sample Python code for retreiving the control input for the \(i\) the state in xgrid:

1p5 = nts_ctrlr[encode3[i], ...]
2p7 = ctlr[p5[2]:p5[2]+p5[0], ...]
3uset = np.argwhere(p7[:, 0] == q).squeeze()
4if(uset.size > 1):
5       uid = rng.choice(uset, 1)  # randomly pick one
6else:
7       uid = int(uset)
8       u = U[p7[uid, 1], :].squeeze()

The DBA state also needs to be update during each sampling time, and this is where you should use q_prime. Let n_dba be the number of states in the DBA and q be the current DBA state. Then the following line of code updates the DBA state:

q = q_prime[p5[1]*n_dba+q]

Some Python APIs are provided in ROCS to convenient simulation such as the ones for

  • loading/reading control synthesis results or specifications: read_controller_itvl_from_h5, read_controller_itvl_from_mat, read_controller_abst_from_h5, and read_spec_from_txt.

  • converting values to indices: index_in_interval_array, index_in_grid.

  • simulation: simulate_itvl_dba_control, simulate_abst_dba_control

You may also use these APIs if they are actually helpful to you.