traffic.c

2车道交通流模拟元胞自动机换道模型代码,是学习换道模型的基础.

// 2 Roads Merging NaSch CA model
// By David Clarke ( [e-mail] )
// Original source available at www.maths.tcd.ie/~dave/traffic/
//Compile with cc -o traffic traffic.c

#include<stdio.h>
#include<stdlib.h>
#include<math.h>

#define leadin_d 500 
#define overlap_d 10
#define leadout_d 500
#define max_length 2100  // greater then leadin+overlap+leadout
#define speed_limit1_d 5
#define speed_limit2_d 5
#define steps_d 20000 		//Number of Steps updates
#define cars1_d 0		// cars to start with (0 is empty road)
#define cars2_d 0
#define cars0_d 0 // leadout cars
#define prob_d 0.5		// probability of stopping
#define slow_start_d 0.0
#define relax_d 600		// steps ignore before measuring
#define flux_time_d 2000		//steps flux measured over
#define print_d 0
#define split 10000000
#define once 0 			// if 1, run for only one density, for scatter plot.
#define max_density 1 	// limiter 0->1, can stop program for a lower density 

typedef struct {
  int car;
  int vel;
  int moved;
  int pre_vel;			//previous velocity, used in slow to start
  int jtime;
  int origin;
  int num;
} road_t;
FILE *out1;
FILE *outScatter1, *outScatter2;
FILE *out_jtime1, *out_jtime2;

//Function definitions
void print_road (int);
void print_jams ();
void print_flow_density();
void print_scatter (int c, int tstep);
void setup (int);
void setup_leadout ();
void setspeed (int);
void update (int c, int tstep);
void measure_density();
void variables (int argc, char *argv[]);

//global variables
road_t Road1[max_length],Road2[max_length];
int leadin,overlap,leadout, speed_limit[3], steps, cars[3], relax, density_width, flux_time,tot_length,print;
double prob, slow_start;
int count[3],density_start[3],density_stop[3], time_count, density[3], jtimeMax[3],jtimeMin[3];

int main (int argc, char *argv[])
{		/* Start up */
  int i,j, n, d,pr;
  road_t *roadP;
  variables (argc, argv);	//take in variables from command line
  tot_length=leadin+overlap+leadout;
//density measurement parameters.
  density_start[1]=density_start[2]=leadin/2;
  density_start[0]=leadin+overlap+leadout/2;
  density_stop[1]=density_stop[2]=leadin;
  density_stop[0]=leadin+overlap+leadout;
  out1 = fopen ("density_flow.dat", "w");
  outScatter1 = fopen ("trafficYScatter1.dat", "w");
  outScatter2 = fopen ("trafficYScatter2.dat", "w");
  out_jtime1 = fopen ("jtime1.dat", "w");
  out_jtime2 = fopen ("jtime2.dat", "w");
  srand (time (0));		// random seed


while(cars[1]<leadin*max_density&&cars[2]<leadin*max_density){ // Looping for increasing density 
		/* Set Up */
	n = 0;
	setup (1);			// initial positioning of cars
	setup (2);			// initial positioning of cars
	setup_leadout();
	count[1]=count[2]=1;
	if(print){	print_road(1);       // initial print
		print_road(2);       // initial print
	}
		/* Relax */	
	while(n < relax) {n++; //allow system to settle
		setspeed(1); 
		update(1,-1); 
		setspeed(2); 
		update(2,-1);}   

			//	Give each car a number:
	for(i=0;i<leadin+overlap+leadout;i++){
		if(Road1[i].car==1){Road1[i].num=count[Road1[i].origin];
			count[Road1[i].origin]++;}
		if(Road2[i].car==1){Road2[i].num=count[Road2[i].origin];
			count[Road2[i].origin]++;}
		}

	time_count = flux_time;
	count[0]=count[1]=count[2] = n = 0;
	density[0]=density[1]=density[2]=0;
  		
		/* Run Model */
  while (n < steps) {
	n++;
	time_count--;
	setspeed (1);   
	update (1,n); 
	setspeed (2);   
	update (2,n);
	if(print){printf("r1"); print_road(1);
		 	printf("r2"); print_road(2);
	}

	if(once){
			print_scatter(1,n); // Print Scatter diagram
			print_scatter(2,n); 
	}
	measure_density();
	if (time_count == 0)print_flow_density();
	
	//    print_jams ();
  } // steps
  if(once) break; // do program just once (for scatter)
	cars[1]+=5;
	cars[2]+=5;
	cars[0]+=5;
}
  fclose (out1);
  fclose (outScatter1);
  fclose (outScatter2);
  fclose (out_jtime1);
  fclose (out_jtime2);
fprintf(stderr,"\nprogram complete\n");
exit (0);

}		//end of main


void setup (int c){
	int i;
	int cars_used;
	road_t *roadP;
	switch(c){
		case 1: roadP=&Road1[0];
		break;
		case 2:  roadP=&Road2[0];
		break;
		default: fprintf(stderr,"Error! Unknown case in setup\n");}
		cars_used=cars[c];
	for (i = 0; i < tot_length; i++) {
		(roadP+i)->car = 0;
		(roadP+i)->vel = 0; }

	while (cars_used > 0) {
    	for (i = 0; i < leadin; i++) {	// initialise array
      		if (((int) (0 + (1.0*leadin/cars_used)*(double)rand()
				/(RAND_MAX + 1.0)))==1
				&& cars_used > 0 && (roadP+i)->car == 0) {
			(roadP+i)->car = 1;
			(roadP+i)->vel =(int)(0+(speed_limit[c]+1)*(double)rand()/(RAND_MAX + 1.0));
			cars_used--;
			(roadP+i)->jtime=-split;
			(roadP+i)->origin=c;
      		}
	}}
}
void setup_leadout (){
	int i;
	int cars_used=cars[0];
	road_t *roadP;
	roadP=&Road1[leadin+overlap];

	while (cars_used > 0) {
    	for (i = 0; i < leadout; i++) {	// initialise array
      		if (((int) (0 + (1.0*leadout/cars_used)*(double)rand()
				/(RAND_MAX + 1.0)))==1
				&& cars_used > 0 && (roadP+i)->car == 0) {
			(roadP+i)->car = 1;
			(roadP+i)->vel =(int)(0+(speed_limit[1]+1)*(double)rand()/(RAND_MAX + 1.0));
			cars_used--;
			(roadP+i)->jtime=-split;
			if(cars[1]==0)
			(roadP+i)->origin=2;
			else if(cars[2]==0)
			(roadP+i)->origin=1;
			else if((int)((1+(1.0*cars[2])/cars[1])*(double)rand()/(RAND_MAX +1.0))==0)
			(roadP+i)->origin=1;
			else
			(roadP+i)->origin=2;

//fprintf(stderr,"rand: %d\n",(int)((1+(1.0*cars[2])/cars[1])*(double)rand()/(RAND_MAX +1.0)));
      		}
	}}
}
void
setspeed (int c){
	int i, j, dist, k, d, p,length;
	road_t *roadP, *roadO;
	switch(c){
		case 1: length=leadin+overlap+leadout;
		roadP=&Road1[0];
		break;
		case 2: length=leadin+overlap;
		roadP=&Road2[0];
		break;
		default: fprintf(stderr,"Error! Unknown case in setspeed\n");
	}

	for (i = 0; i < length; i++) {	// for each road segment
    	if ((roadP+i)->car == 1) {	// only look at segments with cars
			(roadP+i)->pre_vel = (roadP+i)->vel;
			//1. accelerate
			if ((roadP+i)->vel < speed_limit[c])(roadP+i)->vel++; 
			
			dist = 0;
			// Find distance to next car (up to speed limit)
			for(j=1;j<=speed_limit[c];j++){
				if (i+j>=length) {
					if(c==1){  //Main branch -off the edge
					// Find average 
						if((roadP+i)->origin==1)roadO=&Road1[0];
						else roadO=&Road2[0];
						if((roadO+i+j-length)->car==0) dist++;
						else { 	break;}
					}
					else  break;   // Slip lane, stop at end
				}			//loop road
				else if ((roadP+i+j)->car == 0) dist++;
				else break;
 			}

			//debug
			//printf("i:%d j:%d k:%d dist:%d\n",i,j,k,dist);
			//printf("vel:%d\n",(roadP+i)->vel);

			//2. Deceleration (with infinate breaking power)
			if ((roadP+i)->vel > dist) (roadP+i)->vel = dist;				
			//printf("vel:%d\n",Road[i].vel);
			
			// 3. Random Deceleration
      		if((roadP+i)->vel>0&&((int)(0+(1.0/prob)*(double)rand()/
				(RAND_MAX+1.0)))==0) 
				(roadP+i)->vel--;		
		}
	}
}

void update (int c,int tstep){
	int i, n = 0, d,length,v,k,m;
	road_t *roadP, *roadO;
	switch(c){
		case 1: length=leadin+overlap+leadout;
		roadP=&Road1[0];
		break;
		case 2: length=leadin+overlap;
		roadP=&Road2[0];
		break;
		default: fprintf(stderr,"Error! Unknown case in setspeed\n");
	}
	roadO=roadP;
	for(i=0;i<length;i++)(roadP+i)->moved=0;// Don't move cars twice
	
	for(i=length-1;i>=0;i--){ //move backwards because of merging lane
		if ((roadP+i)->car==1&&(roadP+i)->moved==0){ //if car&not moved
			k=m=0;

			/* ### Merging  ### */
			/* ###  Rules   ### */

/* Move over at the earliest opportunity, then move up lane as much as possible ie. min(vel,dist to next car) */

			if(c==2&&i>=leadin){	//in merging zone
				v=Road2[i].vel;
				for(k=0;k<=v;k++){  //what k can we move over at?
					if(Road1[i+k].car==0&&Road1[i+k-1].car==0){
						m++; //succeeded in moving over (and 1 space to spare)
						d=k; } // max advancement
					if(m>0&&Road1[i+k].car==1) break; //moved over and behind next car
				}
				if(m==0){ //couldn't move, move in lane at reduced speed
					if(v>1){
						v--;
						if(Road2[i+v].car==1)fprintf(stderr,"Crash while not merging!!");
						Road2[i+v].car=1;
						Road2[i].car=0;
						Road2[i+v].vel=v;
						Road2[i+v].jtime=Road2[i].jtime+1;
						Road2[i+v].origin=2;
						Road2[i+v].num=Road2[i].num;
						}
				}
				else{
					if(Road1[i+d].car==1)fprintf(stderr,"Crash while merging!!");
					Road1[i+d].car=1;
					Road1[i+d].vel=v;
					Road1[i+d].jtime=Road2[i].jtime+1;
					Road1[i+d].num=Road2[i].num;
					Road1[i+d].origin=2;
					Road2[i].car=0;
				}
			}
	/* ### Other motion ### */
	
			else if((roadP+i)->vel == 0)(roadP+i)->jtime++; // car doesn't move
			else if(i+(roadP+i)->vel>=length){ // car loops.
				count[0]++; 
				d=(roadP+i)->vel+i-length;//car now loops around
				if((roadP+i)->origin==2){ roadO=&Road2[0];}
					/* calculate its journey time */
					// if time <0 car came from middle,ignore 
				if((roadP+i)->jtime>0&&tstep!=-1){
					if((roadP+i)->origin==1)
						fprintf(out_jtime1,"1\t%d\t%d\n",tstep,(roadP+i)->jtime);
					else fprintf(out_jtime2,"2\t%d\t%d\n",tstep,(roadP+i)->jtime);
				}
				if((roadO+d)->car == 1)fprintf(stderr,"\n *** ### Car crash! ### ***\n");
				(roadO+d)->car = 1;
				(roadO+d)->vel = (roadP+i)->vel;
				(roadO+d)->jtime = 0;
				(roadO+d)->num = (roadP+i)->num;
				(roadO+d)->origin = (roadP+i)->origin;
				(roadO+d)->moved = 1;
				(roadP+i)->car = 0;
			}
				// 4. Update & Move
			else {
				d = (roadP+i)->vel + i;
				
				if((roadP+d)->car == 1)fprintf(stderr,"\n *** ### Car crash! ### ***\n");
				(roadP+d)->car = 1;
				(roadP+d)->vel = (roadP+i)->vel;
				(roadP+d)->jtime = (roadP+i)->jtime+1;
				(roadP+d)->num = (roadP+i)->num;
				(roadP+d)->origin = (roadP+i)->origin;
				(roadP+d)->moved = 1;
				(roadP+i)->car = 0;
				if(i<leadin&&d>=leadin){count[c]++;}//car has left the zone
			} 
			
		}
	}
}  //end of update()


void measure_density(){
	int i,j;
	road_t *roadP;
	for(j=0;j<3;j++){
		switch(j){
			case 0: roadP=&Road1[0]; break;
			case 1: roadP=&Road1[0]; break;
			case 2: roadP=&Road2[0]; break;
			default: fprintf(stderr,"error measuring density");}
	for (i=density_start[j];i<density_stop[j];i++){
		if ((roadP+i)->car == 1) density[j]++;
      	}  	
//	if(print)  printf("D%d:%d ",j,density[j]);
	}
}
void
print_road (int c){
	int i,length;
	road_t *roadP;
	switch(c){
		case 1: roadP=&Road1[0];
		length=leadin+overlap+leadout;
		break;
		case 2:  roadP=&Road2[0];
  		length=leadin+overlap;
  		break;
		default: fprintf(stderr,"Error! Unknown case in print road\n");}

	for (i = 0; i < tot_length; i++) {
		if(i==0||i==leadin||i==leadin+overlap||i==leadin+overlap+leadout)
	  	printf("|");
    	if ((roadP+i)->car == 0) {
			printf ("-");
    	}
	    else {
      printf ("%d", (roadP+i)->vel);    // print velocities
//      printf ("%d", (roadP+i)->origin); // print origin
//      printf ("%d", (roadP+i)->num%10); //print car number mod 10
    	}
 	}
	printf ("\n");
}

void print_scatter(int c,int tstep){
	int i,length;
	road_t *roadP;
	FILE *Scatter;

	switch(c){
		case 1: roadP=&Road1[0];
		length=leadin+overlap+leadout;
		Scatter=outScatter1;
		break;
	case 2:  roadP=&Road2[0];
		length=leadin+overlap;
		Scatter=outScatter2;
		break;
		default:fprintf(stderr,"Error! Unknown case in print scatter\n");}

	for (i = 0; i < length; i++) {
    	if ((roadP+i)->car == 1&&(roadP+i)->num%5==0) {
			if(c==2)fprintf(Scatter,"%d\t%d\n",tstep,i);
	/* Separate for plotting different colours based on origin */
		else if(Road1[i].origin==1)fprintf(Scatter,"%d\t%d\t-\n",tstep,i);
		else if(Road1[i].origin==2)fprintf(Scatter,"%d\t-\t%d\n",tstep,i);
		else fprintf(stderr,"Error, origin not 1 or 2");
		}
	}
}

void print_flow_density(){
int j;
	for(j=0;j<3;j++){
		fprintf(out1,"%lf\t%lf\t",(density[j]*1.0)/(flux_time*(density_stop[j]-density_start[j])),(count[j]*1.0)/flux_time);
		time_count = flux_time;
		count[j] = 0;
		density[j] = 0;
	}
	fprintf(out1,"\n");
}

//Take in variables from command line
void variables (int argc, char *argv[]){
	if (argc == 1) {
		leadin = leadin_d;
		overlap = overlap_d;
		leadout = leadout_d;
		speed_limit[1] = speed_limit1_d;
		speed_limit[2] = speed_limit2_d;
		steps = steps_d;
		cars[1] = cars1_d;
		cars[2] = cars2_d;
		cars[0] = cars0_d;
		relax = relax_d;
		flux_time = flux_time_d;
		prob = prob_d;
		slow_start = slow_start_d;
		leadin = leadin_d;
		print= print_d;
	}
	else if(argc !=14 ) {
		fprintf (stderr,"\nUsage: %s <leadin> <overlap> <leadout> <cars1> <cars2> <cars leadout> <speed1> <speed2> <steps> <relax> <prob> <fluxTime> <print>\n\t \n\n",argv[0]);
		fprintf(stderr,"Try: %s  %d %d %d %d %d %d %d %d %d %d %lf %d %d\n\n", argv[0],leadin_d,overlap_d,leadout_d,cars1_d,cars2_d,cars0_d,speed_limit1_d,speed_limit2_d,steps_d,relax_d,prob_d,flux_time_d,print_d);
		exit (1);
    }
	else {
	if (sscanf (argv[1], "%d", &leadin) != 1) {
		fprintf(stderr,"\narg1: road length not an int\n\n");
		exit (1);}
	if (sscanf (argv[2], "%d", &overlap) != 1) {
		fprintf(stderr,"\narg2: road length not an int\n\n");
		exit (1); }
	if (sscanf (argv[3], "%d", &leadout) != 1) {
		fprintf(stderr,"\narg3: road length not an int\n\n");
		exit (1); }
	if(leadin+overlap+leadout>max_length){
		fprintf(stderr,"\n Total length greater then max length\n\n");
		exit(1);}
	if (sscanf (argv[4], "%d", &cars[1]) != 1) {
		fprintf(stderr,"\narg4:Cars1 not an int\n\n");
		exit (1); }
	if (sscanf (argv[5], "%d", &cars[2]) != 1) {
		fprintf(stderr,"\narg5:Cars2 not an int\n\n");
		exit (1); }
	if (sscanf (argv[6], "%d", &cars[0]) != 1) {
		fprintf(stderr,"\narg5:Cars leadout not an int\n\n");
		exit (1);}
	if (sscanf (argv[7], "%d", &speed_limit[1]) != 1) {
		fprintf(stderr,"\narg7:Speed Limit 1 not an int\n\n");
		exit (1); }
	if (sscanf (argv[8], "%d", &speed_limit[2]) != 1) {
		fprintf(stderr,"\narg8:Speed Limit 2 not an int\n\n");
		exit (1);  }
	if (sscanf (argv[9], "%d", &steps) != 1) {
		fprintf(stderr,"\narg9: Steps not an int\n\n");
		exit (1); }
	if (sscanf (argv[10], "%d", &relax) != 1) {
		fprintf(stderr,"\narg10: relax not an int\n\n");
		exit (1);}
	if (sscanf (argv[11], "%lf", &prob) != 1 || prob > 1) {
		fprintf(stderr,"\narg11:Slowing prob not less then 1\n\n");
		exit (1);}
	if(cars[1]>leadin||cars[2]>leadin){
		fprintf(stderr,"\nToo many cars for road length\n\n");
		exit(1);}
	if (sscanf (argv[12], "%d", &flux_time) != 1) {
		fprintf(stderr,"\narg12:Flux time not an int\n\n");
		exit (1);}
	if (sscanf (argv[13], "%d", &print) != 1&&(print==0||print==1)) {
		fprintf(stderr,"\narg13:Print is not 0 or 1\n\n");
		exit (1);}
    }
}