Below is the explaination regarding this code from Jordan. Attached is the code in its last form.
Hi Lucas,
I'll try to fit a bunch of stuff into this email, so hopefully some of it is useful/interesting.
As far as the evolutionary algorithm goes, I'll send you the code and let you look at it yourself. the paper that I am developing it from comes from https://www.researchgate.net/publication/282857432_INTEGRATION_OF_REMCOM_XFDTD_AND_OCTAVE_FOR_THE_AUTOMATED_OPTIMAL_DESIGN_OF_MICROWAVE_ANTENNAS. My code has a couple differences in that I didn't look at the past k iterations like the paper does, but instead I have 5 parallel designs being run and look at how many of those improve the output. This is subtly different and might not be good, because it considers each design at the same time, so the mean doesn't have time to adjust before it is reevaluated if that makes sense. Anyways, just something to think about.
The basic structure for the code is that it is run with command line arguments, so that you compile it with whatever name ( I usually call it evolved_dipole, but doesn't matter). So it is run as
$./evolved_dipole --start
for the first time in a run and then every subsequent time in a run
$./evolve_dipole --cont
The --start will create the initial population, and record the parameters in a .csv file called handshake.csv. The --cont will read in a file called handshook.csv that theoretically will have the output of the antenna simulations for each antenna.
The first obvious thing I can think of that is missing in this script is that it doesn't write a number to a txt file in the watch folder, but I'll explain that later. The second obvious thing that I didn't add to this is the checking of the exit condition d/do < some value (see paper if this is confusing). The third thing I can think of is that I don't have any constraints on the values that the script produces. It will probably be valuable to include some constraints on these values at the very least so that you don't get negative distances. In addition, this script should be easily generalizable by increasing NVAR and changing the mean and deviation vectors. The code is also not particularly pretty so I apologize for that. I tried to comment a decent amount.
Then, there is the XF script. So the XF script should input the antenna parameters, create the simulations and then start them. Then output the data. One of the things that I never ended up doing here is that the data output can only be done after the simulations have finished running, so you'll need to figure out how to get that to work if you use XF, so I'll include the scripts separately. For the output.xmacro script you will need to specify the simulation and run id each time it is used. I think it might be possible to just have a while function that will wait x seconds and then reevaluate and the condition is whether the data is available, but that might not be the best way.
Then, we get to the controlling portion of the code. I have a batch script which should control a bash script (batch submission for clusters is weird). so the bash script theoretically should load the xf module, run evolved_dipole --start, watch a folder and see if a certain number is written to a file, run the xf simulation, watch again for a different number, run the evolved_dipole --cont (or arasim first eventually) and then it creates a loop and should run theoretically until the exit condition is reached in which case --cont will need to write a different value than before and everything should end (i don't know that I've included this part yet).
The big problem here is that calling the XF script from the command line is more difficult than it originally appeared to be. According to Remcom (the company that creates XF), you should be able to run an XF script from the command line with the option --execute-macro-script=/path/to/macro. However, this isn't supported (we think) by the version of XF that we have on OSC, and so they are looking in (talk to Carl abou this) to updating XF and how much that would cost. I'm not entirely sure that this is a solution either, because this requires calling the GUI for XF and I'm not sure that this is able to be done in a batch submission (clusters don't like using cluster computing with GUIs). Thus, it might be worthwhile to look into using NEC or HFSS, which I don't know anything about.
Have fun and let me know if you need any clarification or help,
Jordan |
//Macro to build a dipole antenna made of PEC in XF
//J. Potter -- 06/17
var path = "/users/PAS0654/osu8742/Evolved_Dipole/handshake.csv";
var file = new File(path);
file.open(1);
var hands = file.readAll();
//hands=hands.replace(/,/g, '');
//hands=hands.replace(/[(]/g, '');
//hands=hands.replace(/[)]/g, '');
//Output.println(hands);
var radii=[];
var lengths=[];
//Only take the lines with functions
var lines = hands.split('\n');
for(var i = 0;i < lines.length;i++){
if(i>=8)
{
var params = lines[i].split(",");
radii[i-8]=params[0];
lengths[i-8]=params[1];
Output.println(radii[i-8]);
Output.println(lengths[i-8]);
}
}
file.close();
// set variables for length of the dipoles, radius of the dipoles and units for the lengths
var dipoleLength = 10;
var dipoleRadius = 5;
var connectorLength=1.0;
var gridSize = .25;
var units = " cm" //don't delete the space
var frequency = .5;
//App.saveCurrentProjectAs("C:/Users/Jordan/Desktop/Emacs_Scripts/XFScripts/XFDipoleTest1.xf");
function CreateDipole()
{
App.getActiveProject().getGeometryAssembly().clear();
//create a new sketch
var dipole = new Sketch();
var base = new Ellipse( new Cartesian3D(0,0,0), new Cartesian3D( dipoleRadius + units,0,0 ), 1.0, 0.0, Math.PI*2 );
dipole.addEdge(base);
//add depth to the circle
var extrude = new Extrude( dipole, dipoleLength + units );
//create a recipe and model -- still need to figure what a recipe is...
var dipoleRecipe = new Recipe();
dipoleRecipe.append(extrude);
var dipoleModel = new Model();
dipoleModel.setRecipe(dipoleRecipe);
dipoleModel.name = "Dipole Antenna Test";
//set locations of the left and right dipole segments
dipoleModel.getCoordinateSystem().translate( new Cartesian3D(0,0,1 + units));
var dipoleInProject1 = App.getActiveProject().getGeometryAssembly().append(dipoleModel);
dipoleModel.getCoordinateSystem().translate(new Cartesian3D(0,0,(-1 - dipoleLength) + units));
var dipoleInProject2 = App.getActiveProject().getGeometryAssembly().append(dipoleModel);
// Now set the material for the Dipole:
var pecMaterial = App.getActiveProject().getMaterialList().getMaterial( "PEC" );
if( null == pecMaterial )
{
Output.println( "\"PEC\" material was not found, could not associate with the antenna." );
}
else
{
App.getActiveProject().setMaterial( dipoleInProject1, pecMaterial );
App.getActiveProject().setMaterial( dipoleInProject2, pecMaterial );
}
//zoom to view the extent of the creation
View.zoomToExtents();
}
function CreatePEC() //borrowed from XF demo
{
//Make the material. We will use PEC, or Perfect Electrical Conductor:
var pec = new Material();
pec.name = "PEC";
var pecProperties = new PEC(); // This is the electric properties that defines PEC
var pecMagneticFreespace = new MagneticFreespace(); // We could make a material that acts as PEC and PMC, but in this case we just care about electrical components.
var pecPhysicalMaterial = new PhysicalMaterial();
pecPhysicalMaterial.setElectricProperties( pecProperties );
pecPhysicalMaterial.setMagneticProperties( pecMagneticFreespace );
pec.setDetails( pecPhysicalMaterial );
// PEC is historically a "white" material, so we can easily change its appearance:
var pecBodyAppearance = pec.getAppearance();
var pecFaceAppearance = pecBodyAppearance.getFaceAppearance(); // The "face" appearance is the color/style associated with the surface of geometry objects
pecFaceAppearance.setColor( new Color( 255, 255, 255, 255 ) ); // Set the surface color to white. (255 is the maximum intensity, these are in order R,G,B,A).
// Check for an existing material
if( null != App.getActiveProject().getMaterialList().getMaterial( pec.name ) )
{
App.getActiveProject().getMaterialList().removeMaterial( pec.name );
}
App.getActiveProject().getMaterialList().addMaterial( pec );
}
function CreateAntennaSource()
{
// Here we will create our waveform, create our circuit component definition for the feed, and create
// a CircuitComponent that will attach those to our current geometry.
var waveformList = App.getActiveProject().getWaveformList();
// clear the waveform list
waveformList.clear();
var parameterList = App.getActiveProject().getParameterList();
parameterList.clear();
parameterList.addParameter("freq",".5");
// Create a sinusoidal input wave
var waveform = new Waveform();
var sine = new RampedSinusoidWaveformShape();
//sine.setFrequency("freq" + " GHz"); //Uncomment if using parameter sweep
sine.setFrequency(frequency + " GHz"); //Comment if no parameter sweep
waveform.setWaveformShape( sine );
waveform.name ="Sinusoid";
var waveformInList = waveformList.addWaveform( waveform );
// Now to create the circuit component definition:
var componentDefinitionList = App.getActiveProject().getCircuitComponentDefinitionList();
// clear the list
componentDefinitionList.clear();
// Create our Feed
var feed = new Feed();
feed.feedType = Feed.Voltage; // Set its type enumeration to be Voltage.
// Define a 50-ohm resistance for this feed
var rlc = new RLCSpecification();
rlc.setResistance( "50 ohm" );
rlc.setCapacitance( "0" );
rlc.setInductance( "0" );
feed.setImpedanceSpecification( rlc );
feed.setWaveform( waveformInList ); // Make sure to use the reference that was returned by the list, or query the list directly
feed.name = "50-Ohm Voltage Source";
var feedInList = componentDefinitionList.addCircuitComponentDefinition( feed );
// Now create a circuit component that will be the feed point for our simulation
var componentList = App.getActiveProject().getCircuitComponentList();
componentList.clear();
var component = new CircuitComponent();
component.name = "Source";
component.setAsPort( true );
// Define the endpoints of this feed - these are defined in world position, but you can also attach them to edges, faces, etc.
var coordinate1 = new CoordinateSystemPosition( 0, 0,0 );
var coordinate2 = new CoordinateSystemPosition( 0, 0, "1" + units );
component.setCircuitComponentDefinition( feedInList );
component.setEndpoint1( coordinate1 );
component.setEndpoint2( coordinate2 );
componentList.addCircuitComponent( component );
}
function CreateGrid()
{
// Set up the grid spacing for the dipole antenna
var grid = App.getActiveProject().getGrid();
var cellSizes = grid.getCellSizesSpecification();
cellSizes.setTargetSizes( Cartesian3D( gridSize + units, gridSize + units, gridSize + units ) );
// And we need to set the Minimum Sizes - these are the minimum deltas that we will allow in this project.
// We'll use the scalar ratio of 20% here.
cellSizes.setMinimumSizes( Cartesian3D( ".5", ".5", ".5" ) );
cellSizes.setMinimumIsRatioX( true );
cellSizes.setMinimumIsRatioY( true );
cellSizes.setMinimumIsRatioZ( true );
grid.specifyPaddingExtent( Cartesian3D( "20", "20", "20" ), Cartesian3D( "20", "20", "20" ), true, true );
}
function CreateSensors()
{
// Here we will create a sensor definition and attach it to a near-field sensor on the surface of one of our objects.
var sensorDataDefinitionList = App.getActiveProject().getSensorDataDefinitionList();
sensorDataDefinitionList.clear();
// Create a sensor
var farSensor = new FarZoneSensor();
farSensor.retrieveTransientData = true;
farSensor.setAngle1IncrementRadians(Math.PI/12.0);
farSensor.setAngle2IncrementRadians(Math.PI/12.0);
farSensor.name = "Far Zone Sensor";
var FarZoneSensorList = App.getActiveProject().getFarZoneSensorList();
FarZoneSensorList.clear();
FarZoneSensorList.addFarZoneSensor( farSensor );
}
function CreateAntennaSimulationData()
{
// This function modifies the NewSimulationData parameters of the project.
// They're called "New" because they get applied every time we create an instance of this
// project and place it on our simulation queue.
var simData = App.getActiveProject().getNewSimulationData();
// These should already be set, however just to make sure the current project is set up correctly
simData.excitationType = NewSimulationData.DiscreteSources;
simData.enableSParameters = false;
// Set convergence threshold to -40 dB, and maximum number of timesteps to 10000.
var terminationCriteria = simData.getTerminationCriteria();
terminationCriteria.convergenceThreshold = -40;
terminationCriteria.setMaximumSimulationTime( "20000*timestep" );
// Do not attempt to collect steady-state data
simData.getFOIParameters().collectSteadyStateData = true;
// Construct parameter sweep
//var sweep = simData.getParameterSweep();
//sweep.parameterName = "freq";
//sweep.setParameterValuesByCount(.1, 2, 5); //(start value ,end value, # of steps)
//simData.enableParameterSweep = true;
}
function QueueSimulation()
{
// Create and start the simulation. Project must be saved in order to run.
var simulation = App.getActiveProject().createSimulation( true );
Output.println( "Successfully created the simualtion." );
var projectId = simulation.getProjectId();
var simulationId = simulation.getSimulationId();
var numRuns = simulation.getRunCount();
}
for(var i = 0;i < 5;i++){
var dipoleLength = lengths[i];
var dipoleRadius = radii[i];
//Actually create the material and then the dipole
CreatePEC();
CreateDipole();
CreateAntennaSource();
CreateGrid();
CreateSensors();
CreateAntennaSimulationData();
QueueSimulation();
}
|
/*
* File: Evolved_Dipole.cpp
* Author: Jordan Potter
*
* Created on July 14, 2017, 11:07 AM
*/
#include <cstdlib>
#include <iostream>
#include <cmath>
#include <fstream>
#include <string>
#include <sstream>
using namespace std;
const int NVAR=2; //number of variables (r and l)
const int NPOP=5; //size of population
double rand_normal(double mean, double stddev)
{//Box muller method to generate normal numbers
static double n2 = 0.0;
static int n2_cached = 0;
if (!n2_cached)
{
double x, y, r;
do
{
x = 2.0*rand()/RAND_MAX - 1;
y = 2.0*rand()/RAND_MAX - 1;
r = x*x + y*y;
}
while (r == 0.0 || r > 1.0);
{
double d = sqrt(-2.0*log(r)/r);
double n1 = x*d;
n2 = y*d;
double result = n1*stddev + mean;
n2_cached = 1;
return result;
}
}
else
{
n2_cached = 0;
return n2*stddev + mean;
}
}
void GenerateSample(double new_val[][NVAR], double mean_val[], double std_dev[],int n){
double u;
for(int i=0;i<NVAR;i++)
{
//Create a sample of mean + stdev*(gaussian sample) to use
u=rand_normal(0,.5);
new_val[n][i]=mean_val[i]+std_dev[i]*u;
//Include conditionals for bounds and constraints here
int errors=0;
if(new_val[n][i] <= 0)
{
i--; //If a parameter is less than 0, decrement so this value is created again
errors++;
}
if(errors >= 50) //If there are 50 errors something is going wrong, so stop
{
cout << "Error: Sample Generation Failed" << endl;
return;
}
}
}
void Simulation(double x[][NVAR], double y[NVAR], double mean[], double deviation[], double best[]){
//Write to CSV file. Might need to slightly adjust for generations after #1...
//For now, this will only write the last generation to file
ofstream handshake;
handshake.open("handshake.csv");
handshake << "C++ ESTRA -- Jordan Potter" << "\n";
handshake << "Mean Vector:" << "\n"; //Write the Current Mean Values to handshake.csv
for(int i=0;i<NVAR;i++)
{
if(i==(NVAR-1))
{
handshake << mean[i] << "\n";
}
else
{
handshake << mean[i] << ",";
}
}
handshake << "Deviation Vector:" << "\n"; //Write the Current Deviation Values to handshake.csv
for(int i=0;i<NVAR;i++)
{
if(i==(NVAR-1))
{
handshake << deviation[i] << "\n";
}
else
{
handshake << deviation[i] << ",";
}
}
handshake << "Best Antenna Vector:" << "\n"; //Write the Current Deviation Values to handshake.csv
for(int i=0;i<(NVAR+1);i++)
{
if(i==(NVAR))
{
handshake << best[i] << "\n";
}
else
{
handshake << best[i] << ",";
}
}
handshake << "Generation:" << "\n"; //Write the Generation Values to handshake.csv
for(int i=0;i<NPOP;i++)
{
for(int j=0;j<NVAR;j++)
{
if(j==(NVAR-1))
{
handshake << x[i][j] << "\n";
}
else
{
handshake << x[i][j] << ",";
}
}
}
handshake.close();
//Set Value in another file, so that controller.sh knows to start XF
}
void Read(double x[][NVAR], double y[NVAR], double mean[], double deviations[], double best[]){
//Import values from Handshake for mean/x population/current deviations
ifstream handshake;
handshake.open("handshake.csv");
string csvContent[NPOP+8]; //contain each line of csv
string strToDbl; //gets overwritten to contain each cell of a line. Then transferred to x,y,mean,deviations,best
for(int i=0;i<(NPOP+8);i++)
{
getline(handshake,csvContent[i]); //read in mean
if(i==2)
{
istringstream stream(csvContent[i]);
for(int j=0;j<NVAR;j++)
{
getline(stream,strToDbl,',');
mean[j] = atof(strToDbl.c_str());
}
}
else if(i==4) //read in deviations
{
istringstream stream(csvContent[i]);
for(int j=0;j<NVAR;j++)
{
getline(stream,strToDbl,',');
deviations[j] = atof(strToDbl.c_str());
}
}
else if(i==6) //read in the current best values
{
istringstream stream(csvContent[i]);
for(int j=0;j<NVAR+1;j++)
{
getline(stream,strToDbl,',');
best[j] = atof(strToDbl.c_str());
}
}
else if(i>=8) //read in the generation
{
istringstream stream(csvContent[i]);
for(int j=0;j<NVAR;j++)
{
getline(stream,strToDbl,',');
x[i-8][j] = atof(strToDbl.c_str());
}
}
}
//Import Values from Simulation from Handshook that XF (or eventually ARASim) will write
//May need to adjust if format of handshook changes
ifstream handshook;
handshook.open("handshook.csv");
string strToDbl2[NPOP+2];
for(int i=0;i<(NPOP+2);i++)
{
getline(handshook,strToDbl2[i]);
if(i>=2)
{
y[i-2] = atof(strToDbl2[i].c_str());
}
}
handshook.close();
}
int Mutate(double x[][NVAR],double out[], double best[], double mean[]){
//check to see if the output from x[i][] is larger than the output from x[m].
//i.e. output[i]>mean_output
int count=0;
for(int i=0;i<NPOP;i++)
{
if(out[i]>best[0])
{
best[0]=out[i];
for(int j=0;j<NVAR;j++)
{
best[j+1]=x[i][j];
mean[j]=x[i][j];
}
count++;
}
}
return count;
}
void CheckConvergence(int numSuccess, double p, double q, double d[]){
//If things aren't getting better than expand search radius.
//However, if things are getting better tighten search radius.
//Check this over k? values of x.
//P(output)>p then d=d/q
//if not then d=d*q
double P = (double) numSuccess/NPOP;
if(P >= p)
{
for(int i=0;i<NVAR;i++)
{
d[i]=d[i]/q;
}
}
else
{
for(int i=0;i<NVAR;i++)
{
d[i]=d[i]*q;
}
}
}
int main(int argc, char** argv) {
double m[NVAR]={.5,.5}; //mean values of r and l
double x[NPOP][NVAR]; //produced value of r and l
double d[NVAR] = {.25,.25}; //standard deviation vector
double best[NVAR+1]={-40,0,0}; //array to store param values and output score of top antenna {fit_score,r,l}
double output[NPOP]; //array to store scores of each generations
double q = .9; //factor to adjust search radius by
double p = .2; //ratio of samples that must be correct to adjust search radius
int numSuccess; //number of samples that improve on the best value
srand(time(0));
if(argc != 2)
cout << "Error: Usage. Specify --start or --cont" << endl;
else
{
if(string(argv[1]) == "--start")
{
for(int i=0;i<NPOP;i++)
{
GenerateSample(x,m,d,i); //Generate NPOP samples
}
Simulation(x, output, m, d, best); // Write current population to disk as well as mean,deviations and best population
}
else if(string(argv[1]) == "--cont")
{
//Read Function (read in population to x[],mean to m[],deviations to d[], best to best[] and scores to output[].
Read(x,output,m,d,best);
numSuccess = Mutate(x, output, best, m);
CheckConvergence(numSuccess, p, q, d);
//Check if d/d_o<lambda has been achieved. Consider adding this to check convergence function.
cout << "The best parameter values are:" << endl;
for(int i =0; i <NVAR; i++)
{
cout << m[i] << endl; //Consider writing to a log file
}
for(int i=0;i<NPOP;i++)
{
GenerateSample(x,m,d,i);
}
Simulation(x, output, m, d, best);
}
else
{
cout << "Error: Specify --start or --cont" << endl;
}
}
return 0;
}
|