Commit f7e78f5c authored by ga85rup's avatar ga85rup

modifications for continuous input added

parent 731d69fa
File added
File mode changed from 100644 to 100755
......@@ -14,7 +14,6 @@ using namespace std;
const int sDIM = 10; /* System dimension */
const int iDIM = 2; /* Input dimension */
const double T = 25; /* Sampling time */
size_t P = 3; /* Cardinality of (quantized) input set */
size_t N = 12; /* Temporal Horizon */
typedef std::array<double,sDIM> state_type;
......@@ -24,19 +23,9 @@ template<class F>
void ode_solver(F rhs, state_type &x, size_t nint, double h);
/* ODE */
auto system_post = [](state_type &x, int u) -> void {
/* Mapping abstract input to system inputs */
double us[iDIM]={0};
if(u==0)
{us[0]=0;us[1]=0;}
if(u==1)
{us[0]=0;us[1]=1;}
if(u==2)
{us[0]=1;us[1]=0;}
auto system_post = [](state_type &x, double* u) -> void {
/* ode describing thermal model*/
auto rhs=[us](state_type &xx, const state_type &x) -> void {
auto rhs=[u](state_type &xx, const state_type &x) -> void {
const double a=0.05;
const double ae2=0.005;
const double ae5=0.005;
......@@ -45,10 +34,10 @@ auto system_post = [](state_type &x, int u) -> void {
const double te=12;
const double th=100;
xx[0]=(-a-ae)*x[0]+a*x[1]+ae*te;
xx[1]=(-4*a-ae2-ah*us[0])*x[1]+a*x[0]+a*x[6]+a*x[8]+a*x[2]+ae2*te+ah*th*us[0];
xx[1]=(-4*a-ae2-ah*u[0])*x[1]+a*x[0]+a*x[6]+a*x[8]+a*x[2]+ae2*te+ah*th*u[0];
xx[2]=(-2*a-ae)*x[2]+a*x[1]+a*x[3]+ae*te;
xx[3]=(-2*a-ae)*x[3]+a*x[2]+a*x[4]+ae*te;
xx[4]=(-4*a-ae5-ah*us[1])*x[4]+a*x[3]+a*x[7]+a*x[5]+a*x[9]+ae5*te+ah*th*us[1];
xx[4]=(-4*a-ae5-ah*u[1])*x[4]+a*x[3]+a*x[7]+a*x[5]+a*x[9]+ae5*te+ah*th*u[1];
xx[5]=(-a-ae)*x[5]+a*x[4]+ae*te;
xx[6]=(-a-ae)*x[6]+a*x[1]+ae*te;
xx[7]=(-a-ae)*x[7]+a*x[4]+ae*te;
......@@ -56,14 +45,14 @@ auto system_post = [](state_type &x, int u) -> void {
xx[9]=(-a-ae)*x[9]+a*x[4]+ae*te;
};
size_t nint = 5; /* no. of time step for ode solving */
double h=T/nint; /* time step for ode solving (T is an sampling time) */
double h=T/nint; /* time step for ode solving (T is the sampling time) */
ode_solver(rhs,x,nint,h); /* Runga Kutte solver */
};
/* defining constrains for the controller
ul[i] : upper limit for the temperature in ith dimension
ll[i] : lower limit for the temperature in ith dimension */
auto constrain = [](state_type y) -> bool {
/* defining safe set for the controller
ul[i] : upper bound on the temperature in ith room
ll[i] : lower bound on the temperature in ith room */
auto setBounds = [](state_type y) -> bool {
double ul=21.8;
double ll=18;
bool s = true;
......@@ -90,15 +79,15 @@ int main(){
state_type xs;
xs[0]=17;xs[1]=17;xs[2]=17;xs[3]=17;xs[4]=17;xs[5]=17;xs[6]=17;xs[7]=17;xs[8]=17;xs[9]=17;
/* initial state of the system */
state_type x0;
x0[0]=19;x0[1]=19;x0[2]=19;x0[3]=19.8;x0[4]=20.8;x0[5]=19.8;x0[6]=19;x0[7]=19.8;x0[8]=19;x0[9]=19.8;
/* approximation on abstraction */
double epsilon = 0.1;
/* number of time samples the controller should run */
size_t t = 20;
/************************************************************
input information
*************************************************************/
const size_t P = 3; /* Number of elements in input set*/
double ud[P][iDIM]={{0,0},{0,1},{1,0}};
double *UD[P];
for(size_t i=0;i<P;i++){
UD[i]=ud[i];
}
/******************************************************
Symbolic model construction
......@@ -116,7 +105,7 @@ int main(){
/* Computing the constrain set */
std::cout<<"Computing the constrain set ... ";
tt.tic();
BDD set = ab.getConstrainSet(system_post,constrain,xs);
BDD set = ab.getAbstractSet(system_post,setBounds,xs,iDIM,P,UD);
tt.toc();
if(set == ddmgr.bddZero()){
std::cout << "Set is empty !!";
......@@ -140,14 +129,26 @@ int main(){
return 0;
}
/******************************************************
closed-loop simulation
******************************************************/
/* initial state of the system */
state_type x0;
x0[0]=19;x0[1]=19;x0[2]=19;x0[3]=19.8;x0[4]=20.8;x0[5]=19.8;x0[6]=19;x0[7]=19.8;x0[8]=19;x0[9]=19.8;
/* approximation on abstraction */
double epsilon = 0.1;
/* number of time samples the controller should run */
size_t t = 20;
/* find out mode sequence which is very close to the output */
std::cout<<"Computed initial state is...\n";
/* finding inital stae in abstraction for safety specification */
/* finding initial state in abstraction for safety specification */
BDD w0 = ab.getAbsState(system_post,set,x0,xs,epsilon,sDIM);
/* open-loop simulation */
ab.openLoopSim(C,w0,system_post,x0,sDIM,t);
ab.closedLoopSim(C,w0,system_post,x0,sDIM,t);
return 0;
......
File mode changed from 100644 to 100755
File mode changed from 100644 to 100755
......@@ -11,11 +11,10 @@
using namespace std;
const int sDIM = 6; /* System dimension */
const int iDIM = 2; /* Input dimension */
const size_t sDIM = 6; /* System dimension */
const size_t iDIM = 2; /* Input dimension */
const double T = 25; /* Sampling time */
size_t P = 3; /* Cardinality of (quantized) input set */
size_t N = 12; /* Temporal Horizon */
const size_t N = 6; /* Temporal Horizon */
typedef std::array<double,sDIM> state_type;
......@@ -24,19 +23,10 @@ template<class F>
void ode_solver(F rhs, state_type &x, size_t nint, double h);
/* ODE */
auto system_post = [](state_type &x, int u) -> void {
/* Defining system inputs to system*/
double ub[iDIM]={0};
if(u==0)
{ub[0]=0;ub[1]=0;}
if(u==1)
{ub[0]=0;ub[1]=1;}
if(u==2)
{ub[0]=1;ub[1]=0;}
/* ode describing thermal model*/
auto rhs=[ub](state_type &xx, const state_type &x) -> void {
auto system_post = [](state_type &x, double* u) -> void {
/* ode describing six-room thermal model*/
auto rhs=[u](state_type &xx, const state_type &x) -> void {
const double a=0.05;
const double ae1=0.005;
const double ae4=0.005;
......@@ -44,26 +34,26 @@ auto system_post = [](state_type &x, int u) -> void {
const double ah=0.0036;
const double te=12;
const double th=100;
xx[0] = (-3*a-ae1-ah*ub[0])*x[0]+a*x[1]+a*x[2]+a*x[4]+ae1*te+ah*th*ub[0];
xx[0] = (-3*a-ae1-ah*u[0])*x[0]+a*x[1]+a*x[2]+a*x[4]+ae1*te+ah*th*u[0];
xx[1] = (-2*a-ae)*x[1]+a*x[0]+a*x[3]+ae*te;
xx[2] = (-2*a-ae)*x[2]+a*x[0]+a*x[3]+ae*te;
xx[3] = (-3*a-ae4-ah*ub[1])*x[3]+a*x[1]+a*x[2]+a*x[5]+ae4*te+ah*th*ub[1];
xx[3] = (-3*a-ae4-ah*u[1])*x[3]+a*x[1]+a*x[2]+a*x[5]+ae4*te+ah*th*u[1];
xx[4] = (-a-ae)*x[4]+a*x[0]+ae*te;
xx[5] = (-a-ae)*x[5]+a*x[3]+ae*te;
};
size_t nint = 5; /* no. of time step for ode_solving */
double h=T/nint;
ode_solver(rhs,x,nint,h);
double h=T/nint; /* time step for ode solving (T is an sampling time) */
ode_solver(rhs,x,nint,h); /* Runga Kutte solver */
};
/* defining constrains for the controller
ul[i] : upper limit of the ith dimension
ll[i] : lower limit of the ith dimension */
auto constrain = [](state_type y) -> bool {
/* defining target set for the controller
ul[i] : upper bound on target temperature in ith room
ll[i] : upper bound on target temperature in ith room */
auto setBounds = [](state_type y) -> bool {
double ul=21.0;
double ll=17.5;
bool s = true;
for(int j = 0; j < sDIM; j++)
for(size_t j = 0; j < sDIM; j++)
{
if( y[j] >= ul || y[j] <= ll )
{
......@@ -85,25 +75,25 @@ int main(){
/* source state of the system */
state_type xs;
xs[0]=17;xs[1]=17;xs[2]=17;xs[3]=17;xs[4]=17;xs[5]=17;
/* initial state of the system */
state_type x0;
x0[0]=16;x0[1]=16;x0[2]=16;x0[3]=16;x0[4]=16;x0[5]=16;
/* approximation on abstraction */
double epsilon = 0.2;
/* number of time samples the controller should run */
size_t t = 40;
/************************************************************
input information
***********************************************************/
/* lower bounds on inputs */
double lb[sDIM]={0,0};
/* upper bounds on inputs */
double ub[sDIM]={1,1};
/* quantization parameter */
double eta[sDIM]={.5,1};
/******************************************************
Symbolic model construction
******************************************************/
/* defining Symbolic Set for the transition relation */
SymbolicSetSpace ss(ddmgr,P,N);
SymbolicSetSpace ss(ddmgr,iDIM,lb,ub,eta,N);
tt.tic();
ss.addAbsStates();
/* defining abstraction class */
getAbstraction<state_type> ab(&ss);
tt.toc();
......@@ -112,7 +102,7 @@ int main(){
/* Computing the constrain set */
std::cout<<"Computing the constrain set ... ";
tt.tic();
BDD set = ab.getConstrainSet(system_post,constrain,xs);
BDD set = ab.getAbstractSet(system_post,setBounds,xs,iDIM,lb,ub,eta);
tt.toc();
if(set == ddmgr.bddZero()){
......@@ -137,15 +127,28 @@ int main(){
std::cout << "No controller is found !!";
return 0;
}
/******************************************************
closed-loop simulation
******************************************************/
/* initial state of the system */
state_type x0;
x0[0]=16;x0[1]=16;x0[2]=16;x0[3]=16;x0[4]=16;x0[5]=16;
/* approximation on abstraction */
double epsilon = 0.8;
/* number of time samples the controller should run */
size_t t = 40;
/* find out mode sequence which is very close to the output */
std::cout<<"Computed the initial state is...\n";
/* finding inital stae in abstraction for safety specification */
/* finding initial state in abstraction for safety specification */
BDD w0 = ab.getAbsState(system_post,x0,xs,epsilon,sDIM);
/* open-loop simulation */
ab.openLoopSim(C,w0,system_post,x0,sDIM,t);
/* closed-loop simulation */
ab.closedLoopSim(C,w0,system_post,x0,sDIM,t);
return 0;
......
......@@ -8,16 +8,40 @@ ah=3.6e-3;
te=12;
th=100;
Dt=5;
lb=[0,0]
ub=[1,1]
eta=[0.5,1]
%variable Declaration
Y=zeros(sDIM,6);
X1=zeros(sDIM,1);
%control inputs
filename = 'controller.txt';
input=importdata(filename);
for i=1:length(input)
u(1,i)=bitand(floor(input(i)/2),1);
u(2,i)=bitand(input(i),1);
for i=1:length(lb)
nInput(i)=floor(ub(i)/eta(i))-ceil(lb(i)/eta(i))+1;
end
for i=1:length(lb)-1
temp1=1;
for j=i:length(lb)-1
temp1=temp1*nInput(j+1);
end
xx(i)=temp1;
end
xx;
for j=1:length(input)
for i=1:length(lb)
if i<length(lb)
temp2=floor(input(j)/xx(i));
temp3=mod(temp2,xx(i));
u(i,j)=lb(i)+temp3*eta(i);
elseif i==length(lb)
temp4=mod(input(j),xx(i-1));
u(i,j)=lb(i)+temp4*eta(i);
end
end
end
% initial conditions
X_init=[16 16 16 16 16 16];
%X_init=19.5*ones(sDIM,D_tau);
......@@ -52,15 +76,15 @@ plot(Time,Y);
figure(2);
plot(Y(1,:),Y(2,:));
hold on;
plot([17.5, 17.5 ,21,21,17.5],[17.5,21,21,17.5,17.5]);
plot([17.5, 17.5 ,21.5,21.5,17.5],[17.5,21.5,21.5,17.5,17.5]);
figure(3);
plot(Y(3,:),Y(4,:));
hold on;
plot([17.5, 17.5 ,21,21,17.5],[17.5,21,21,17.5,17.5]);
plot([17.5, 17.5 ,21.5,21.5,17.5],[17.5,21.5,21.5,17.5,17.5]);
figure(4);
plot(Y(5,:),Y(6,:));
hold on;
plot([17.5, 17.5 ,21,21,17.5],[17.5,21,21,17.5,17.5]);
plot([17.5, 17.5 ,21.5,21.5,17.5],[17.5,21.5,21.5,17.5,17.5]);
figure(5);
plot(Time1,Uem(1,:));
figure(6);
......
File mode changed from 100644 to 100755
......@@ -19,50 +19,25 @@ size_t N = 12; /* Temporal Horizon */
typedef std::array<double,sDIM> state_type;
/* ODE */
auto system_post = [](state_type &x, int p) -> void {
auto system_post = [](state_type &x, double* u) -> void {
state_type z = x;
/* the ode describing the system */
const double v = 70 ;
const double l = 0.25 ;
const double T = 0.002777;
const double q = 0.25;
double b[iDIM];
if(p==0)
{
b[0] = 6;
b[1] = 0;
b[2] = 8;
b[3] = 0;
b[4] = 0;
}
if(p==1)
{
b[0] = 6;
b[1] = 0;
b[2] = 0;
b[3] = 0;
b[4] = 0;
}
if(p==2)
{
b[0] = 0;
b[1] = 0;
b[2] = 8;
b[3] = 0;
b[4] = 0;
}
x[0] = (1-T*v/l)*z[0] + b[0];
x[1] = T*v/l*z[0] + (1-T*v/l-q) *z[1] + b[1];
x[2] = T*v/l*z[1] + (1-T*v/l) *z[2] + b[2];
x[3] = T*v/l*z[2] + (1-T*v/l) *z[3] + b[3];
x[4] = T*v/l*z[3] + (1-T*v/l) *z[4] + b[4];
x[0] = (1-T*v/l)*z[0] + u[0];
x[1] = T*v/l*z[0] + (1-T*v/l-q) *z[1] + u[1];
x[2] = T*v/l*z[1] + (1-T*v/l) *z[2] + u[2];
x[3] = T*v/l*z[2] + (1-T*v/l) *z[3] + u[3];
x[4] = T*v/l*z[3] + (1-T*v/l) *z[4] + u[4];
};
/* defining constrains for the controller
/* defining safety bounds for the controller
ul[i] : upper limit of the ith dimension
ll[i] : lower limit of the ith dimension */
auto constrain = [](state_type y) -> bool {
auto setBounds = [](state_type y) -> bool {
state_type ul;
state_type ll;
for(int i = 0; i < sDIM; i++)
......@@ -87,11 +62,9 @@ int main(){
/*to measure time */
TicToc tt;
/* CUDD Manager */
Cudd ddmgr;
/* source state of the system */
state_type xs;
xs[0] = 3.8570;
......@@ -99,17 +72,16 @@ int main(){
xs[2] = 3.3750;
xs[3] = 8.5177;
xs[4] = 8.5177;
/* initial state of the system */
state_type x0;
x0[0] = 1.41783;
x0[1] = 4.92788;
x0[2] = 10.6711;
x0[3] = 9.53216;
x0[4] = 14.5308;
double epsilon = 0.00085;
/* number of times the controller should rum */
size_t t = 14;
/************************************************************
input information
*************************************************************/
const size_t P = 3; /* Number of elements in input set*/
double ud[P][iDIM]={{6,0,8,0,0},{6,0,0,0,0},{0,0,8,0,0}};
double *UD[P];
for(size_t i=0;i<P;i++){
UD[i]=ud[i];
}
/******************************************************
Symbolic model construction
......@@ -127,7 +99,7 @@ int main(){
/* Computing the constrain set */
std::cout<<"Computing the constrain set ... ";
tt.tic();
BDD set = ab.getConstrainSet(system_post,constrain,xs);
BDD set = ab.getAbstractSet(system_post,setBounds,xs,iDIM,P,UD);
tt.toc();
if(set == ddmgr.bddZero()){
......@@ -150,12 +122,26 @@ int main(){
return 0;
}
/******************************************************
* closed-loop simulation
******************************************************/
/* initial state of the system */
state_type x0;
x0[0] = 1.41783;
x0[1] = 4.92788;
x0[2] = 10.6711;
x0[3] = 9.53216;
x0[4] = 14.5308;
/* approximation on abstraction */
double epsilon = 0.00085;
/* number of times the controller should run */
size_t t = 14;
std::cout<<"Computing the initial state ... ";
BDD w0 = ab.getAbsState(system_post,set,x0,xs,epsilon,sDIM);
/* to get output */
ab.openLoopSim(C,w0,system_post,xs,sDIM,t);
ab.closedLoopSim(C,w0,system_post,xs,sDIM,t);
return 0;
......
File mode changed from 100644 to 100755
File mode changed from 100644 to 100755
File mode changed from 100644 to 100755
No preview for this file type
File mode changed from 100644 to 100755
File mode changed from 100644 to 100755
File mode changed from 100644 to 100755
......@@ -28,6 +28,8 @@ friend class abstractionMode;
Cudd* mgr_;
size_t P_;
size_t N_;
size_t dim_;
double* eta_;
size_t *nofBddVars_;
size_t **indBddVars_;
size_t *nofAbsStates_;
......@@ -78,6 +80,60 @@ public :
}
}
/* new constructor for continuous input */
SymbolicSetSpace( Cudd &mgr, const size_t dim, const double* lb, const double* ub, const double* eta, const size_t N){
for (size_t i=0; i<dim; i++) {
if((ub[i]-lb[i])<eta[i]) {
std::ostringstream os;
os << "Error: scots::SymbolicSet: each interval must contain at least one cell.";
throw std::invalid_argument(os.str().c_str());
}
}
dim_=dim;
eta_ = new double[dim];
double Nl, Nu;
size_t nInput[dim];
size_t temp1=1;
for(size_t i=0;i<dim;i++){
eta_[i]=eta[i];
Nl=std::ceil(lb[i]/eta[i]);
Nu=std::floor(ub[i]/eta[i]);
nInput[i]=(size_t)std::abs(Nu-Nl)+1;
temp1=temp1*nInput[i];
}
const size_t P=temp1;
mgr_ = &mgr;
P_ = P;
N_ = 2*N + 1;
nvars_ = 0;
nofBddVars_ = new size_t [N_];
indBddVars_ = new size_t* [N_];
nofAbsStates_ = new size_t [N_];
symbolicSetSpace_ = mgr.bddZero();
for(size_t j = 0; j < N_; j++)
{
nofAbsStates_[j] = P;
if( P == 1 )
nofBddVars_[j] = 1;
else
nofBddVars_[j] = std::ceil(log2(P));
indBddVars_[j] = new size_t [nofBddVars_[j]];
for(size_t i = 0; i < nofBddVars_[j]; i++)
{
BDD var = mgr.bddVar();
indBddVars_[j][i] = var.NodeReadIndex();
}
}
for(size_t i = 0; i < N_; i++)
{
for(size_t j = 0; j < nofBddVars_[i]; j++)
{
nvars_++;
}
}
}
~SymbolicSetSpace(){
delete[] nofAbsStates_;
......
File mode changed from 100644 to 100755
File mode changed from 100644 to 100755
File mode changed from 100644 to 100755
This diff is collapsed.
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment