- Do hundreds of simulations using the conical solver to control the rocket.
- Record the state of the rocket and what the conical solver said at each timestep.
- Train a neural network to predict what the conical solver is going to mandate.

- There are no guarantees that it's going to actually work. It might crash.
- It's not going to be optimal. It might run out of fuel ... and then crash.

Off.

Use keys W,A,D to control the rocket manually. 'c' toggles AINote that this isn't optimal. You can kind-of see that it sort-of works, but it's not perfect at all. It sometimes bounces on the ground on the way to a successful landing.

For an example, here's the entire 3-vector library that I use in my simulation. There are very few languages that could express these functions as clearly and concisely. Even if you would normally use a library for 3-vectors, the point still stands that Kotlin is extremely clear and concise:

class V3(var x:Double=0.0, var y:Double=0.0, var z:Double=0.0){ operator fun minus(a : V3) = V3(x - a.x, y - a.y,z - a.z) operator fun plus(a : V3) = V3(x + a.x, y + a.y,z + a.z) operator fun times(dt: Double) = V3(x*dt,y*dt,z*dt) operator fun rem(a:V3) = Math.sqrt((x-a.x)*(x-a.x) + (y-a.y)*(y-a.y) + (z-a.z)*(z-a.z)) operator fun times(a:V3) = V3(y * a.z - z * a.y, z * a.x - x * a.z, x * a.y - y * a.x) fun length2() = x*x + y*y + z*z fun length() = Math.sqrt(x*x+y*y+z*z) infix fun dot(d: V3) = x*d.x + y * d.y + z*d.z override fun toString() = "($x,$y,$z)" fun copy() = V3(x,y,z) fun addAndScaleIP(d: V3, dt: Double) { x+=d.x * dt y+=d.y * dt z+=d.z * dt } fun addIP(d: V3) { x+=d.x y+=d.y z+=d.z } fun scale(dt: Double) = V3(x*dt,y*dt,z*dt) fun isEqual(v3: V3) = x==v3.x&&y==v3.y&&z==v3.z fun normalise() = scale(1.0 / this.length()) fun isFinite() = x.isFinite() && y.isFinite() && z.isFinite() fun toPov() = "< $x,$y,$z >" fun toDoubleArray() = doubleArrayOf(x,y,z) } fun main(args:Array< String >){ val a = V3(1.2,2.5,-0.2) val b = V3(1.6,-2.5,0.2) println("a=$a") println("b=$b") val c = a + b * 2.0 // scaling val d = a dot c // infix operators val e = a * b + V3(d,0.0,0.0) // cross product println("a + b = $c") println("a + b = ${a+b}") println("e = $e") }

import tensorflow as tf import numpy as np import json import matplotlib.pyplot as plt files = [ "C:/data/tfSimple6dofH.json", "C:/data/tfSimple6dofI.json", "C:/data/tfSimple6dofJ.json", "C:/data/tfSimple6dofK.json", "C:/data/tfSimple6dofL.json", "C:/data/tfSimple6dofM.json", "C:/data/tfSimple6dofN.json", "C:/data/tfSimple6dofO.json", "C:/data/tfSimple6dofP.json" ] datasets = [json.load(open(_)) for _ in files] dataset = [] for ds in datasets: for d2 in ds: dataset.append(d2) print("length of dataset = " + str(len(dataset)) + " rows") cutoff = (len(dataset)*3)//4 train = dataset[:cutoff] test = dataset[cutoff:] x_train = np.array([_['input'] for _ in train]) y_train = np.array([_['output'] for _ in train])*0.01+0.5 x_test = np.array([_['input'] for _ in test]) y_test= np.array([_['output'] for _ in test])*0.01+0.5 model = tf.keras.models.Sequential([ tf.keras.layers.Dense(20, activation=tf.nn.sigmoid), tf.keras.layers.Dense(20, activation=tf.nn.sigmoid), tf.keras.layers.Dense(len(train[0]['output']), activation=tf.nn.sigmoid) ]) model.compile(optimizer='adam', loss='mse'#'sparse_categorical_crossentropy', #metrics=['accuracy'] ) scoreIn = [] scoreOut = [] for i in range(1000): model.fit(x_train, y_train, epochs=5, verbose=10, steps_per_epoch=5) a = model.evaluate(x_train, y_train) b = model.evaluate(x_test, y_test) print("iteration " + str(i) + " score = " + str(int(a*1e5)) + " vs " + str(int(b*1e5))+"\n") scoreIn.append(a) scoreOut.append(b) w=model.get_weights() j = [] for i in range(len(w)//2): j.append({ 'weights': w[i*2+0].tolist(), 'bias': w[i*2+1].tolist() }) json.dump({'network':j}, open("C:/tmp/weights6dofB.json",'w')) plt.plot(scoreOut[10:]);plt.plot(scoreIn[10:])This isn't really complicated at all. It may have been better, in an ideal world, to use Kotlin throughout or python throughout.

I have to say that the documentation, at least in the current version, for ECOS, is hard to use. The C interface is only described in one paragraph that doesn't tie the variable names to the equivalent equations properly, and which overrun the edge of the PDF page.

The way to learn to use it is to start off with the simplest possible problem, and then keep altering it until you get a sensible result, and then keep adding to it until you end up with a list of what's what.

class Problem{ var nVars : Int = 0 val G = mutableListOf< DoubleArray>() val q = mutableListOf< Int>() val c = mutableListOf< Double>() val h = mutableListOf< Double>() var pEquals = 0 var mConeVars = 0 var lOrthants = 0 var sCones = 0 fun save(fn : String){ val ret = mutableListOf< String>() ret.add("(n $nVars)") ret.add("(p $pEquals)") ret.add("(m $mConeVars)") ret.add("(l $lOrthants)") ret.add("(s $sCones)") val gMatrix = newMatrix(G) val gNonZero = gMatrix.nonZeros() ret.add("(g $gNonZero)") ret.add("(end)") ret.add("") for(i in 0 until q.size){ ret.add("(q $i,${q[i]})") } for(i in 0 until c.size){ ret.add("(c $i,${c[i]})") } for(i in 0 until h.size){ ret.add("(h $i,${h[i]})") } for(i in 0 until G.size){ for(j in 0 until G[i].size){ val g2 = G[i][j] if (g2!=0.0){ ret.add("(g $i,$j,$g2)") } } } ret.save(fn) } fun addVariable(cVal : Double):Int{ this.nVars++ this.c.add(cVal) if (G.size>0) throw Exception("Add variables before adding G contributions") return this.nVars-1 } fun addOrthant(g:DoubleArray, hVal:Double){ if (sCones>0) throw Exception("Add orthants before cones") lOrthants++ mConeVars++ G.add(g) h.add(hVal) } fun addCone(g:List< DoubleArray>,hVals:List< Double>){ sCones++ q.add(g.size) mConeVars+=g.size h.addAll(hVals) G.addAll(g) } fun constrainTwice(i: Int, x: Double, tolerance:Double=0.0) { val d = DoubleArray(nVars) d[i] = 1.0 addOrthant(d, x+0.001+tolerance) val d2 = DoubleArray(nVars) d2[i] = -1.0 addOrthant(d2, -x+0.001+tolerance) } } class trajectoryInput( var pos:V3, var vel: V3, var g:Double, var dt:Double, var nSteps:Int, var verticalTime:Double, var maxThrust: Double) class ConeOptions(var justDown : Boolean, var endingFall : Boolean) fun coneTrajectory(inp:trajectoryInput, coneOptions: ConeOptions) : Problem{ /* Variables: velocities in each step abs thrust used in each step change in velocity must be less than abs thrust used. minimise total abs thrust abs thrust must be less than threshold. */ var prob = Problem() var vx = mutableListOf< Int>() var vy = mutableListOf< Int>() var vz = mutableListOf< Int>() var thrust = mutableListOf< Int>() for (i in 0 until inp.nSteps) { vx.add(prob.addVariable(0.0)) vy.add(prob.addVariable(0.0)) vz.add(prob.addVariable(0.0)) if (i < inp.nSteps - 1) { thrust.add(prob.addVariable(1.0)) // might be -1.0 } } val dai = { DoubleArray(prob.nVars) } prob.constrainTwice(vx[0], inp.vel.x) prob.constrainTwice(vy[0], inp.vel.y) prob.constrainTwice(vz[0], inp.vel.z) var driftTolerance = 0.3 prob.constrainTwice(vx[inp.nSteps - 1], 0.0, driftTolerance) prob.constrainTwice(vy[inp.nSteps - 1], 0.0, driftTolerance) prob.constrainTwice(vz[inp.nSteps - 1], 0.0, driftTolerance) if (coneOptions.endingFall) { // for the last few steps, I only want vertical movement. val verticalSteps = Math.min(inp.nSteps, (inp.verticalTime / inp.dt).toInt()) val acceleration = inp.maxThrust + inp.g // g negative, so this'll be about 20m/s normally. for (i in inp.nSteps - verticalSteps until inp.nSteps - 1) { prob.constrainTwice(vx[i], 0.0, driftTolerance) prob.constrainTwice(vz[i], 0.0, driftTolerance) if (i == inp.nSteps - verticalSteps) { val timeToEnd = (inp.nSteps - i) * inp.dt val height = -0.5 + 0.5 * timeToEnd * timeToEnd * acceleration val dMinH = dai() dMinH[vy[i]] = -0.5 for (j in i until inp.nSteps - 1) { dMinH[vy[j]] = -1.0 } // dMinH is now the height at this timestep prob.addOrthant(dMinH * -1.0, height * -1.0) } } } // The sum of the velocities must be equal to the negative of initial position, so that it ends up at 0. val dx = dai() val dy = dai() val dz = dai() dx[vx[0]] = inp.dt * 0.5 dy[vy[0]] = inp.dt * 0.5 dz[vz[0]] = inp.dt * 0.5 for (i in 1 until inp.nSteps - 1) { dx[vx[i]] = inp.dt dy[vy[i]] = inp.dt dz[vz[i]] = inp.dt } dx[vx[inp.nSteps - 1]] = inp.dt * 0.5 dy[vy[inp.nSteps - 1]] = inp.dt * 0.5 dz[vz[inp.nSteps - 1]] = inp.dt * 0.5 var posTolerance = 2.0 prob.addOrthant(dx, -inp.pos.x + posTolerance) prob.addOrthant(dy, -inp.pos.y+ posTolerance) prob.addOrthant(dz, -inp.pos.z+ posTolerance) prob.addOrthant(dx * -1.0, inp.pos.x+ posTolerance) prob.addOrthant(dy * -1.0, inp.pos.y+ posTolerance) prob.addOrthant(dz * -1.0, inp.pos.z+ posTolerance) for (i in 0 until inp.nSteps - 1) { val d = dai() d[thrust[i]] = 1.0 prob.addOrthant(d, inp.maxThrust) prob.addOrthant(d*-1.0, 0.0) // this constrains it to be positive, shouldn't be necessary. if (coneOptions.justDown) { // The thrust has to be upward: no boosting downwards. val d2 = dai() d2[vy[i]] = 1.0 d2[vy[i + 1]] = -1.0 prob.addOrthant(d2, -inp.g * inp.dt * 1) } } for (i in 0 until inp.nSteps - 1) { val drad = dai() drad[thrust[i]] = -inp.dt val dx = dai() dx[vx[i]] = 1.0 dx[vx[i + 1]] = -1.0 val dy = dai() dy[vy[i]] = 1.0 dy[vy[i + 1]] = -1.0 val dz = dai() dz[vz[i]] = 1.0 dz[vz[i + 1]] = -1.0 prob.addCone(listOf(drad, dx, dy, dz), listOf(0.0, 0.0, -inp.g*inp.dt, 0.0)) } return prob } class TrajectoryOutput(val pos : List< V3>, val vel : List< V3>, val thrust : List< V3>) fun getTrajUnknownN(inp:trajectoryInput, coneOptions: ConeOptions) : TrajectoryOutput{ var min = 2 var max = 250 var mid = 75 for(i in 0 until 10){ inp.nSteps = mid //println("trying ${inp.nSteps} steps" ) val ret = getTraj(inp, coneOptions) if (ret == null){ min = mid mid = (max+min)/2 if (mid == min) mid++ } else{ max = mid mid = (max + min)/2 if (mid==min || mid==max) return ret } } return getTraj(inp, coneOptions)!! }

class ConicControlOutput(val thrustDir : V3, val pointDir:V3, val expectedPosition : List< V3>){ fun getControl(s : Dynamics.State, d : Dynamics) : DoubleArray{ var pd = pointDir.copy().normalise() var minCosTheta = 0.7 if (pd.y < minCosTheta){ pd.y = 0.0 pd = pd.normalise() * Math.sqrt(1 - minCosTheta*minCosTheta) pd.y=minCosTheta } var target = pd *-1.0 var targetW = (s.rot.fo * target) + s.w*0.3 val ret = DoubleArray(d.thrusters.size) for(i in 0 until d.thrusters.size){ val thruster = d.thrusters[i] val dw = s.rot.objToAbs(thruster.dir * thruster.pos) dot targetW var result = 0.0 if (dw > 0) result = Math.min(dw*0.03,1.0) ret[i] = result } var mainThrust = thrustDir dot (s.rot.objToAbs(d.thrusters[0].dir)) mainThrust *= 1.0 / (d.thrusters[0].force / s.getMass()) if (mainThrust < 0) mainThrust=0.0 if (mainThrust > 1) mainThrust = 1.0 ret[0] = mainThrust return ret } }

## Quantum Mechanics for programmersA quick demo showing how to make simulations of simple quantum mechanics systems in javascript. |

## Mars LanderAn interactive simulation of a Mars Lander trying to land on Mars. The term "interactive" here means that you have to write the software, in Javascript, that will land the probe. |

© Hugo2015. Session @sessionNumber