Commit 88f14dd9 authored by Matthias Rungger's avatar Matthias Rungger
Browse files

The ode solver is now outside the SymbolicModelGrowthBound.hh

parent 8dc1997c
......@@ -3,7 +3,6 @@
#
#CC = g++
CC = clang++
CXXFLAGS = -g -Wall -Wextra -std=c++11
CXXFLAGS = -Wall -Wextra -std=c++11 -O3 -DNDEBUG
#
......@@ -12,33 +11,25 @@ CXXFLAGS = -Wall -Wextra -std=c++11 -O3 -DNDEBUG
SCOTSROOT = ../../..
SCOTSINC = -I$(SCOTSROOT)/src
#
# boost
#
BOOSTINC = -I/opt/local/include
#
# cudd
#
CUDDPATH = ../../../../cudd-2.5.1
CUDDPATH = ../../../../code/ext/cudd-2.5.1
CUDDINC = -I$(CUDDPATH)/include
CUDDLIBS = -lcudd -lobj -ldddmp -lmtr -lepd -lst -lutil
CUDDLPATH = -L$(CUDDPATH)/cudd -L$(CUDDPATH)/util -L$(CUDDPATH)/mtr -L$(CUDDPATH)/st -L$(CUDDPATH)/dddmp -L$(CUDDPATH)/epd -L$(CUDDPATH)/obj
.PHONY: dcdc
TARGET = dcdc
all: $(TARGET)
%.o:%.cc
$(CC) -c $(CXXFLAGS) $(CUDDINC) $(SCOTSINC) $(BOOSTINC) $< -o $@
$(CC) -c $(CXXFLAGS) $(CUDDINC) $(SCOTSINC) $< -o $@
$(TARGET): $(TARGET).o
$(CC) $(CXXFLAGS) -o $(TARGET) $(TARGET).o $(CUDDLPATH) $(CUDDLIBS)
rm $(TARGET).o
clean:
rm $(TARGET)
rm ./$(TARGET) ./$(TARGET).o
......@@ -10,9 +10,8 @@
*
*/
#include <array>
#include <iostream>
#include <boost/array.hpp>
#include <boost/numeric/odeint.hpp>
#include "cuddObj.hh"
......@@ -27,65 +26,90 @@
#define iDIM 1
/* data types for the ode solver */
typedef boost::array<double,2> state_type;
typedef boost::array<double,1> input_type;
/* the ode describing the system */
auto system_ode = []( const state_type &x , const input_type &u, state_type &dxdt, double /*t */) {
const double r0=1.0 ;
const double vs = 1.0 ;
const double rl = 0.05 ;
const double rc = rl / 10 ;
const double xl = 3.0 ;
const double xc = 70.0 ;
const double b[2]={vs/xl, 0};
double a[2][2];
if(u[0]==1) {
a[0][0] = -rl / xl;
a[0][1] = 0;
a[1][0] = 0;
a[1][1] = (-1 / xc) * (1 / (r0 + rc));
} else {
a[0][0] = (-1 / xl) * (rl + ((r0 * rc) / (r0 + rc))) ;
a[0][1] = ((-1 / xl) * (r0 / (r0 + rc))) / 5 ;
a[1][0] = 5 * (r0 / (r0 + rc)) * (1 / xc);
a[1][1] =(-1 / xc) * (1 / (r0 + rc)) ;
}
typedef std::array<double,2> state_type;
typedef std::array<double,1> input_type;
/* forward declaration of the ode solver */
template<class F>
void ode_solver(F rhs, state_type &x, input_type &u, size_t nint, double h);
/* we integrate the dcdc ode by 0.3 sec (the result is stored in x) */
auto dcdc_post = [](state_type &x, input_type &u) -> void {
/* the ode describing the system */
auto system_ode = [](state_type &dxdt, const state_type &x, const input_type &u) -> void {
const double r0=1.0 ;
const double vs = 1.0 ;
const double rl = 0.05 ;
const double rc = rl / 10 ;
const double xl = 3.0 ;
const double xc = 70.0 ;
const double b[2]={vs/xl, 0};
double a[2][2];
if(u[0]==1) {
a[0][0] = -rl / xl;
a[0][1] = 0;
a[1][0] = 0;
a[1][1] = (-1 / xc) * (1 / (r0 + rc));
} else {
a[0][0] = (-1 / xl) * (rl + ((r0 * rc) / (r0 + rc))) ;
a[0][1] = ((-1 / xl) * (r0 / (r0 + rc))) / 5 ;
a[1][0] = 5 * (r0 / (r0 + rc)) * (1 / xc);
a[1][1] =(-1 / xc) * (1 / (r0 + rc)) ;
}
dxdt[0] = a[0][0]*x[0]+a[0][1]*x[1] + b[0];
dxdt[1] = a[1][0]*x[0]+a[1][1]*x[1] + b[1];
};
size_t nint=5; /* number of intermediate step size */
double h=0.1; /* h* nint = sampling time */
ode_solver(system_ode,x,u,nint,h);
};
dxdt[0] = a[0][0]*x[0]+a[0][1]*x[1] + b[0];
dxdt[1] = a[1][0]*x[0]+a[1][1]*x[1] + b[1];
/* computation of the growth bound (the result is stored in r) */
auto radius_post = [](state_type &r, input_type &u) -> void {
/* the ode to determine the radius of the cell which over-approximates the
* attainable set see: http://arxiv.org/abs/1503.03715v1 */
auto growth_bound_ode = [](state_type &drdt, const state_type &r, const input_type &u) {
/* for the dcdc boost converter the growth bound is simply given by the metzler matrix of the system matrices */
const double r0=1.0 ;
const double rl = 0.05 ;
const double rc = rl / 10 ;
const double xl = 3.0 ;
const double xc = 70.0 ;
double a[2][2];
if(u[0]==1) {
a[0][0] = -rl / xl;
a[0][1] = 0;
a[1][0] = 0;
a[1][1] = (-1 / xc) * (1 / (r0 + rc));
} else {
a[0][0] = (-1 / xl) * (rl + ((r0 * rc) / (r0 + rc))) ;
a[0][1] = ((1 / xl) * (r0 / (r0 + rc))) / 5 ;
a[1][0] = 5 * (r0 / (r0 + rc)) * (1 / xc);
a[1][1] =(-1 / xc) * (1 / (r0 + rc)) ;
}
drdt[0] = a[0][0]*r[0]+a[0][1]*r[1];
drdt[1] = a[1][0]*r[0]+a[1][1]*r[1];
};
size_t nint=5; /* number of intermediate step size */
double h=0.1; /* h* nint = sampling time */
ode_solver(growth_bound_ode,r,u,nint,h);
};
/* the ode to determine the radius of the cell which over-approximates the
* attainable set see: http://arxiv.org/abs/1503.03715v1 */
auto growth_bound_ode = []( const state_type &r , const input_type &u, state_type &drdt, double /*t */) {
/* for the dcdc boost converter the growth bound is simply given by the metzler matrix of the system matrices */
const double r0=1.0 ;
const double rl = 0.05 ;
const double rc = rl / 10 ;
const double xl = 3.0 ;
const double xc = 70.0 ;
double a[2][2];
if(u[0]==1) {
a[0][0] = -rl / xl;
a[0][1] = 0;
a[1][0] = 0;
a[1][1] = (-1 / xc) * (1 / (r0 + rc));
} else {
a[0][0] = (-1 / xl) * (rl + ((r0 * rc) / (r0 + rc))) ;
a[0][1] = ((1 / xl) * (r0 / (r0 + rc))) / 5 ;
a[1][0] = 5 * (r0 / (r0 + rc)) * (1 / xc);
a[1][1] =(-1 / xc) * (1 / (r0 + rc)) ;
}
drdt[0] = a[0][0]*r[0]+a[0][1]*r[1];
drdt[1] = a[1][0]*r[0]+a[1][1]*r[1];
};
int main() {
......@@ -137,9 +161,8 @@ int main() {
/****************************************************************************/
scots::SymbolicModelGrowthBound<state_type,input_type> abstraction(&ss, &is);
/* compute the transition relation */
double samplingTime=0.5;
tt.tic();
abstraction.template computeTransitionRelation(system_ode, growth_bound_ode, samplingTime);
abstraction.computeTransitionRelation(dcdc_post, radius_post);
std::cout << std::endl;
tt.toc();
/* get the number of elements in the transition relation */
......@@ -171,4 +194,27 @@ int main() {
return 1;
}
template<class F>
void ode_solver(F rhs, state_type &x, input_type &u, size_t nint, double h) {
/* runge kutte order 4 */
state_type k[4];
state_type tmp;
for(size_t t=0; t<nint; t++) {
rhs(k[0],x, u);
for(size_t i=0;i<sDIM;i++)
tmp[i]=x[i]+h/2*k[0][i];
rhs(k[1],tmp, u);
for(size_t i=0;i<sDIM;i++)
tmp[i]=x[i]+h/2*k[1][i];
rhs(k[2],tmp, u);
for(size_t i=0;i<sDIM;i++)
tmp[i]=x[i]+h*k[2][i];
rhs(k[3],tmp, u);
for(size_t i=0; i<sDIM; i++)
x[i] = x[i] + (h/6)*(k[0][i] + 2*k[1][i] + 2*k[2][i] + k[3][i]);
}
}
......@@ -21,7 +21,7 @@ specification
3. simulate the closed loop in Matlab
a) you need to have the mexfile compile (see the readme in $(SCOTSROOT)/mfiles/mexfiles )
a) you need to compile the mexfile first (see the readme in $(SCOTSROOT)/mfiles/mexfiles )
b) in Matlab add the path and run the simulation
......
......@@ -3,38 +3,29 @@
#
#CC = g++
CC = clang++
CXXFLAGS = -g -Wall -Wextra -std=c++11 -O3
CXXFLAGS = -Wall -Wextra -std=c++11 -O3 -DNDEBUG
#
# scots
#
SCOTSROOT = ../../..
SCOTSINC = -I$(SCOTSROOT)/src
#
# boost
#
BOOSTINC = -I/opt/local/include
#
# cudd
#
CUDDPATH = ../../../../cudd-2.5.1
CUDDPATH = ../../../../code/ext/cudd-2.5.1
CUDDINC = -I$(CUDDPATH)/include
CUDDLIBS = -lcudd -lobj -ldddmp -lmtr -lepd -lst -lutil
CUDDLPATH = -L$(CUDDPATH)/cudd -L$(CUDDPATH)/util -L$(CUDDPATH)/mtr -L$(CUDDPATH)/st -L$(CUDDPATH)/dddmp -L$(CUDDPATH)/epd -L$(CUDDPATH)/obj
TARGET = vehicle
all: $(TARGET)
%.o:%.cc
$(CC) -c $(CXXFLAGS) $(CUDDINC) $(SCOTSINC) $(BOOSTINC) $< -o $@
$(CC) -c $(CXXFLAGS) $(CUDDINC) $(SCOTSINC) $< -o $@
$(TARGET): $(TARGET).o
$(CC) $(CXXFLAGS) -o $(TARGET) $(TARGET).o $(CUDDLPATH) $(CUDDLIBS)
......
A path planning problem for a vehicle
A path planning problem for a vehicle using the bycicle model
1. compile the vehicle.cc file
......@@ -22,12 +22,12 @@ A path planning problem for a vehicle
3. simulate the closed loop in Matlab
a) you need to have the mexfile compile (see the readme in $(SCOTSROOT)/mfiles/mexfiles )
a) you need to compile the mexfile first (see the readme in $(SCOTSROOT)/mfiles/mexfiles )
b) in Matlab add the path and run the simulation
>> addpath(genpath('$(SCOTSROOT)/mfiles'))
>> dcdc
>> vehicle
4. information on the example are found in
......
......@@ -10,9 +10,8 @@
*
*/
#include <array>
#include <iostream>
#include <boost/array.hpp>
#include <boost/numeric/odeint.hpp>
#include "cuddObj.hh"
......@@ -27,38 +26,43 @@
#define sDIM 3
#define iDIM 2
/* functions to setup the state space, input space and obstacles of the vehicle example */
scots::SymbolicSet vehicleCreateStateSpace(Cudd &mgr);
scots::SymbolicSet vehicleCreateInputSpace(Cudd &mgr);
void vehicleCreateObstacles(Cudd &mgr, scots::SymbolicSet &obs);
/* data types for the ode solver */
typedef boost::array<double,3> state_type;
typedef boost::array<double,2> input_type;
/* the ode describing the vehicle */
auto vehicle_ode = []( const state_type &x , const input_type &u, state_type &dxdt, double /*t */) {
double alpha=std::atan(std::tan(u[1])/2.0);
dxdt[0] = u[0]*std::cos(alpha+x[2])/std::cos(alpha);
dxdt[1] = u[0]*std::sin(alpha+x[2])/std::cos(alpha);
dxdt[2] = u[0]*std::tan(u[1]);
typedef std::array<double,3> state_type;
typedef std::array<double,2> input_type;
/* forward declaration of the ode solver */
template<class F>
void ode_solver(F rhs, state_type &x, input_type &u, size_t nint, double h);
/* we integrate the vehicle ode by 0.3 sec (the result is stored in x) */
auto vehicle_post = [](state_type &x, input_type &u) -> void {
/* the ode describing the vehicle */
auto rhs =[](state_type& xx, const state_type &x, input_type &u) -> void {
double alpha=std::atan(std::tan(u[1])/2.0);
xx[0] = u[0]*std::cos(alpha+x[2])/std::cos(alpha);
xx[1] = u[0]*std::sin(alpha+x[2])/std::cos(alpha);
xx[2] = u[0]*std::tan(u[1]);
};
size_t nint=5; /* number of intermediate step size */
double h=0.06; /* h* nint = sampling time */
ode_solver(rhs,x,u,nint,h);
};
/* the ode to determine the radius of the cell which over-approximates the
* attainable set see: http://arxiv.org/abs/1503.03715v1 */
auto growth_bound_ode = []( const state_type &r , const input_type &u, state_type &drdt, double /*t */) {
/* drdt = L r
*
* | 0 0 c |
* L = | 0 0 c |
* | 0 0 0 |
*/
/* computation of the growth bound (the result is stored in r) */
auto radius_post = [](state_type &r, input_type &u) -> void {
double c = std::abs(u[0]*std::sqrt(std::tan(u[1])*std::tan(u[1])/4.0+1));
drdt[0] = c*r[2];
drdt[1] = c*r[2];
drdt[2] = 0;
r[0] = r[0]+c*r[2]*0.3;
r[1] = r[1]+c*r[2]*0.3;
};
/* forward declaration of the functions to setup the state space
* input space and obstacles of the vehicle example */
scots::SymbolicSet vehicleCreateStateSpace(Cudd &mgr);
scots::SymbolicSet vehicleCreateInputSpace(Cudd &mgr);
void vehicleCreateObstacles(Cudd &mgr, scots::SymbolicSet &obs);
int main() {
......@@ -102,16 +106,14 @@ int main() {
/* construct SymbolicSet for the input space */
/****************************************************************************/
scots::SymbolicSet is=vehicleCreateInputSpace(mgr);
std::cout << std::endl << "Input space details:" << std::endl;
/****************************************************************************/
/* setup class for symbolic model computation */
/****************************************************************************/
scots::SymbolicModelGrowthBound<state_type,input_type> abstraction(&ss, &is);
/* compute the transition relation */
double samplingTime=0.3;
tt.tic();
abstraction.template computeTransitionRelation(vehicle_ode, growth_bound_ode, samplingTime);
abstraction.computeTransitionRelation(vehicle_post, radius_post);
std::cout << std::endl;
tt.toc();
/* get the number of elements in the transition relation */
......@@ -236,4 +238,29 @@ scots::SymbolicSet vehicleCreateInputSpace(Cudd &mgr) {
}
template<class F>
void ode_solver(F rhs, state_type &x, input_type &u, size_t nint, double h) {
/* runge kutte order 4 */
state_type k[4];
state_type tmp;
for(size_t t=0; t<nint; t++) {
rhs(k[0],x, u);
for(size_t i=0;i<sDIM;i++)
tmp[i]=x[i]+h/2*k[0][i];
rhs(k[1],tmp, u);
for(size_t i=0;i<sDIM;i++)
tmp[i]=x[i]+h/2*k[1][i];
rhs(k[2],tmp, u);
for(size_t i=0;i<sDIM;i++)
tmp[i]=x[i]+h*k[2][i];
rhs(k[3],tmp, u);
for(size_t i=0; i<sDIM; i++)
x[i] = x[i] + (h/6)*(k[0][i] + 2*k[1][i] + 2*k[2][i] + k[3][i]);
}
}
......@@ -46,7 +46,7 @@ while(1)
u=controller.getInputs(y(end,:));
v=[v; u(1,:)];
[t x]=ode45(@unicycle_ode,[0 .3], y(end,:), odeset('abstol',1e-4,'reltol',1e-4),u(1,:));
[t x]=ode45(@unicycle_ode,[0 .3], y(end,:), [],u(1,:));
y=[y; x(end,:)];
end
......@@ -82,8 +82,6 @@ plot(y(1,1),y(1,1),'.','color',colors(5,:),'markersize',20)
box on
axis([-.5 10.5 -.5 10.5])
%set(gcf,'paperunits','centimeters','paperposition',[0 0 16 10],'papersize',[16 10])
end
......
......@@ -3,38 +3,29 @@
#
#CC = g++
CC = clang++
CXXFLAGS = -g -Wall -Wextra -std=c++11 -O3
CXXFLAGS = -Wall -Wextra -std=c++11 -O3 -DNDEBUG
#
# scots
#
SCOTSROOT = ../../..
SCOTSINC = -I$(SCOTSROOT)/src
#
# boost
#
BOOSTINC = -I/opt/local/include
#
# cudd
#
CUDDPATH = ../../../../cudd-2.5.1
CUDDPATH = ../../../../code/ext/cudd-2.5.1
CUDDINC = -I$(CUDDPATH)/include
CUDDLIBS = -lcudd -lobj -ldddmp -lmtr -lepd -lst -lutil
CUDDLPATH = -L$(CUDDPATH)/cudd -L$(CUDDPATH)/util -L$(CUDDPATH)/mtr -L$(CUDDPATH)/st -L$(CUDDPATH)/dddmp -L$(CUDDPATH)/epd -L$(CUDDPATH)/obj
TARGET = vehicle
all: $(TARGET)
%.o:%.cc
$(CC) -c $(CXXFLAGS) $(CUDDINC) $(SCOTSINC) $(BOOSTINC) $< -o $@
$(CC) -c $(CXXFLAGS) $(CUDDINC) $(SCOTSINC) $< -o $@
$(TARGET): $(TARGET).o
$(CC) $(CXXFLAGS) -o $(TARGET) $(TARGET).o $(CUDDLPATH) $(CUDDLIBS)
......
A path planning problem for a vehicle
A path planning problem for a vehicle using the bycicle model
1. compile the vehicle.cc file
......@@ -22,12 +22,12 @@ A path planning problem for a vehicle
3. simulate the closed loop in Matlab
a) you need to have the mexfile compile (see the readme in $(SCOTSROOT)/mfiles/mexfiles )
a) you need to compile the mexfile first (see the readme in $(SCOTSROOT)/mfiles/mexfiles )
b) in Matlab add the path and run the simulation
>> addpath(genpath('$(SCOTSROOT)/mfiles'))
>> dcdc
>> vehicle
4. information on the example are found in
......
......@@ -10,15 +10,13 @@
*
*/
#include <array>
#include <iostream>
#include <boost/array.hpp>
#include <boost/numeric/odeint.hpp>
#include "cuddObj.hh"
#include "SymbolicSet.hh"
#include "SymbolicModelGrowthBound.hh"
//#include "SymbolicModelsparse.hh"
#include "TicToc.hh"
#include "FixedPoint.hh"
......@@ -28,38 +26,44 @@
#define sDIM 3
#define iDIM 2
/* functions to setup the state space and input space of the vehicle example */
scots::SymbolicSet vehicleCreateStateSpace(Cudd &mgr);
scots::SymbolicSet vehicleCreateInputSpace(Cudd &mgr);
/* data types for the ode solver */
typedef boost::array<double,3> state_type;
typedef boost::array<double,2> input_type;
/* the ode describing the vehicle */
auto vehicle_ode = []( const state_type &x , const input_type &u, state_type &dxdt, double /*t */) {
double alpha=std::atan(std::tan(u[1])/2.0);
dxdt[0] = u[0]*std::cos(alpha+x[2])/std::cos(alpha);
dxdt[1] = u[0]*std::sin(alpha+x[2])/std::cos(alpha);
dxdt[2] = u[0]*std::tan(u[1]);
typedef std::array<double,3> state_type;
typedef std::array<double,2> input_type;
/* forward declaration of the ode solver */
template<class F>
void ode_solver(F rhs, state_type &x, input_type &u, size_t nint, double h);
/* we integrate the vehicle ode by 0.3 sec (the result is stored in x) */
auto vehicle_post = [](state_type &x, input_type &u) -> void {
/* the ode describing the vehicle */
auto rhs =[](state_type& xx, const state_type &x, input_type &u) -> void {
double alpha=std::atan(std::tan(u[1])/2.0);
xx[0] = u[0]*std::cos(alpha+x[2])/std::cos(alpha);
xx[1] = u[0]*std::sin(alpha+x[2])/std::cos(alpha);
xx[2] = u[0]*std::tan(u[1]);
};
size_t nint=5; /* number of intermediate step size */
double h=0.06; /* h* nint = sampling time */
ode_solver(rhs,x,u,nint,h);
};
/* the ode to determine the radius of the cell which over-approximates the
* attainable set see: http://arxiv.org/abs/1503.03715v1 */
auto growth_bound_ode = []( const state_type &r , const input_type &u, state_type &drdt, double /*t */) {
/* drdt = L r
*
* | 0 0 c |
* L = | 0 0 c |
* | 0 0 0 |
*/
/* computation of the growth bound (the result is stored in r) */
auto radius_post = [](state_type &r, input_type &u) -> void {
double c = std::abs(u[0]*std::sqrt(std::tan(u[1])*std::tan(u[1])/4.0+1));
drdt[0] = c*r[2];
drdt[1] = c*r[2];
drdt[2] = 0;