Log-periodic dipole array antenna optimization with SGD

LPDA-SGD.nc

//
// Apply stochastic gradient descent over antenna dimensions array d[]
// to minimize the VSWR.
//
// Some lazy TODO items:
// - implement limits for dimensions in order to prevent model twisting
// - make it more suitable for tiny gaps (requires reviw of the gradient math)
// - stop when the cost function were saturated (or improve the learning rate)
//
// In case of anything, please message me at twitter.com/@mishkathebear
//

real d[50];			// antenna dimensions, big to fit complex models
real dJ[50];			// gradient over the d[] array
real minfreq;
real maxfreq;

real prng_seed;

// Nearly realistic model of a wire wound log-periodic dipole array antenna
model("LPDA") 
{
	real h;			// antenna height
	real wd;		// wire diameter
	int n;			// number of elements
	
	real offset;
	int i, v;
	element feed;

	n = 8;			// as may be derived from other calculators
	h = 0;			// the antenna is in the free space
	wd = 0.0017;		// wire diameter, measured for S = 2.5 mm^2
	
	// params:
	// d[0] - dipole 1 half size (usually the highest freq)
	// d[1] - the boom balanced line wire-to-wire distance
	// d[i] - dipole half length, starting from i = 2
	// d[i+1] - relative offset on the boom
	// d[n] - the closing stub length
	
	feed = wire(0,  0, h-d[1]/2, 0, 0, h+d[1]/2, wd, 3);
	voltageFeed(feed, 1.0, 0.0);
	
	// element 1 (usually of the highest freq)
	wire(0, 0, h-d[1]/2, 0,  d[0], h-d[1]/2, wd, 15);
	wire(0, 0, h+d[1]/2, 0, -d[0], h+d[1]/2, wd, 15);
	
	// generate the rest N-1 elements
	offset = 0;
	i = 2;
	repeat (n - 1) {
		// vertical placement vector
		if ((i / 2) % 2 == 1)
			v = -1;
		else
			v = 1;

		// d[i]   - dipole half length
		// d[i+1] - relative offset on the boom
		wire(offset,        0, h-d[1]/2, offset+d[i+1],     0, h-d[1]/2, wd, 15);
		wire(offset,        0, h+d[1]/2, offset+d[i+1],     0, h+d[1]/2, wd, 15);
		wire(offset+d[i+1], 0, h-v*d[1]/2, offset+d[i+1],  d[i], h-v*d[1]/2, wd, 15);
		wire(offset+d[i+1], 0, h+v*d[1]/2, offset+d[i+1], -d[i], h+v*d[1]/2, wd, 15);
		
		offset = offset + d[i+1];
		
		i = i + 2;
	}
	
	// the closing stub
	wire(offset,      0, h-d[1]/2, offset+d[i], 0, h-d[1]/2, wd, 7);
	wire(offset,      0, h+d[1]/2, offset+d[i], 0, h+d[1]/2, wd, 7);
	wire(offset+d[i], 0, h-d[1]/2, offset+d[i], 0, h+d[1]/2, wd, 3);
	
	freespace();
}

/////////////////////////////////////////////////////////////////////////////

int model_init()
{
	int i;
	int n;
	
	n = 8;				// number of LPDA elements; the same as the
					// n in the model(). Should be global?
	
	minfreq = 2200.0;		// start below 2.4 GHz to take more on it
	maxfreq = 6000.0;		// ... and cover up to 5.8G WiFi

	d[0] = c / maxfreq / 6;		// el. 1, keep it smaller than lambda/4
	d[1] = 0.004;			// usually something about 2*wd to 3*wd
	
	// build the rest N-1 dipoles as linear function
	i = 2;
	repeat(n - 1) {
		// the further on the boom - the bigger
		d[i] = d[i-2] * 1.2;	// something which results to the
					// c/minfreq/4 for the biggest element
		d[i+1] = d[0];		// aligned to high f, impacts directivity
		
		i = i + 2;
	}
	
	d[i] = 0.002;			// length of the closing stub

	return i + 1;			// number of the antenna dimensions
}

/////////////////////////////////////////////////////////////////////////////

int sign(real x)
{
	if (x < 0)
		return -1;
		
	return 1;
}

real prng()
{
	prng_seed = 4001 * prng_seed % 4999;
	return prng_seed / 4999;
}

/////////////////////////////////////////////////////////////////////////////

control()
{
	real J, Jbase, a;
	real f;
	int ndims;
	int epochs;
	int i, k;
	
	prng_seed = 2;		// a carefully picked absolutely random number ;-)
	
	ndims = model_init();

	epochs = 300;		// number of attempts to reach the goal
	a = c/maxfreq / 15;	// the learning rate; try to estimate it automatically
	Jbase = 99.9;		// start from a terrific error and improve it
	
	k = 0;			// the epochs counter
	repeat (epochs) {
		// random pick a frequency from the sweep range, and set on it
		f = minfreq + prng() * (maxfreq - minfreq);
		setFrequency(f);		
		
		// convergence optimization: repeat in the minibatch of two
		// attempts as controlled by the baseline VSWR value (Jbase)
		repeat(2) {
			runModel();
			Jbase = vswr(1);
		
			printf("\nEp. %d: f = %.1fMHz, VSWR=%.2f - ", k, f, Jbase);
		
			// don't optimize the already good result - rather give
			// more time to other freqs
			if (Jbase < 1.5) {
				printf("looks good");
				break;
			}
	
			// calculate the gradient dJ
			i = 0;
			repeat (ndims) {
				// in order to calculate the hyphothezis function h(d[i])
				// we have to temporarily cut the d[i] by small amount of
				// a and evaluate the new model; for big models this will
				// take a while, sorry
				d[i] = d[i] - a * 0.1;
				runModel();
				// recover the model to its original state!
				d[i] = d[i] + a * 0.1;
			
				// the cost function J = vswr(1). Then, log(J) can be used
				// as global component of the gradient (please note, for
				// the ideal case the dJ(vswr=1) = 0), and the log(Jbase/J)
				// is the local part - it shows how much does this particular
				// element affects the cost function.  Also, it is always
				// directed towards minimum, so the gradient too.
				J = vswr(1);
				dJ[i] = log(Jbase) * log(Jbase / J);
			
				printf(" [%d]J=%.2f:dJ=%.4f", i, J, dJ[i]);
				
				i = i + 1;
			}
	
			// now update all dimensions, simultaneously
			i = 0;
			repeat (ndims) {
				d[i] = d[i] - a * dJ[i];
				printf(" d[%d]=%f", i, d[i]);
				i = i + 1;
			}
			
			// one pass done, but do we really need to do it once more?
			if (Jbase < 2) {
				printf(" - was not bad, continue");
				break;
			}
		}

		printf("\n");	

		k = k + 1;
	}

	printf( "\nSGD finished. Evaluating the final model... " ) ;
	frequencySweep(minfreq*0.1, maxfreq*1.1, 51);
	runModel();
	printf("DONE.\n") ;
}