From eaedfaf29c078782b614a2e2467c11c398ec72b2 Mon Sep 17 00:00:00 2001 From: Martin Asprusten Date: Wed, 1 Apr 2026 15:41:25 +0200 Subject: [PATCH] Much guessing --- src/calculations/orbit-calculations.ts | 521 ++++++++++++++++++------- src/gui/worker.ts | 2 +- 2 files changed, 372 insertions(+), 151 deletions(-) diff --git a/src/calculations/orbit-calculations.ts b/src/calculations/orbit-calculations.ts index edd428c..8ff68ba 100644 --- a/src/calculations/orbit-calculations.ts +++ b/src/calculations/orbit-calculations.ts @@ -555,6 +555,97 @@ export function getSpeed(trueAnomaly: number, orbit: Orbit, planet: Body): numbe return Math.sqrt(planet.gravitationalParameter * (1 + 2 * orbit.eccentricity * Math.cos(trueAnomaly) + orbit.eccentricity**2) / orbit.semiLatusRectum); } +export function getOrbitalPeriod(semiMajor: number, body: Body) { + return Math.sqrt(semiMajor**3 / body.gravitationalParameter); +} + +export function perform2dGradientDescent(functionToMinimize: ((variableOne: number, variableTwo: number) => number | null), initialGuessVariableOne: number, initialGuessVaribaleTwo: number, maxDifference?: number, maxTries?: number): [number, number] | null { + let variableOne = initialGuessVariableOne; + let variableTwo = initialGuessVaribaleTwo; + let bestValue = functionToMinimize(variableOne, variableTwo); + + if (bestValue == null) { + return null; + } + + maxDifference = maxDifference ?? 0.5; + maxTries = maxTries ?? 1000; + + // Do some gradient descent on the best values found so far + let tries = 0; + while (bestValue != null && bestValue > maxDifference && tries < maxTries) { + tries++; + + let dfx = functionToMinimize(variableOne+0.0001, variableTwo); + let dfy = functionToMinimize(variableOne, variableTwo+0.0001); + + let estimateDfDx = null; + let estimateDfDy = null; + + if (dfx) { + estimateDfDx = (dfx - bestValue) / 0.0001; + } + + if (dfy) { + estimateDfDy = (functionToMinimize(variableOne, variableTwo+0.0001) ?? bestValue - bestValue) / 0.0001; + } + + let direction: number[][] | null = null; + + if (estimateDfDx && estimateDfDy) { + let slopeNormal = vectorCrossProduct( + [[1], [0], [estimateDfDx]], + [[0], [1], [estimateDfDy]] + ); + + direction = normalizeVector([[slopeNormal[0][0]], [slopeNormal[1][0]], [0]]); + } else if (estimateDfDx) { + direction = [[1], [0], [0]]; + } else if (estimateDfDy) { + direction = [[0], [1], [0]]; + } + + if (!direction) { + break; + } + + let df = functionToMinimize(variableOne + direction[0][0]*0.0001, variableTwo + direction[1][0]*0.0001); + if (!df) { + break; + } + + let estimateDfDv = (df - bestValue) / 0.0001; + + let bestValueImprovement: number = bestValue; + let variableOneUpdate = variableOne; + let variableTwoUpdate = variableTwo; + let divisor = 0; + let deadEnd = false; + + while (bestValueImprovement >= bestValue) { + variableOneUpdate = variableOne - bestValue * direction[0][0] / (estimateDfDv * 2**divisor); + variableTwoUpdate = variableTwo - bestValue * direction[1][0] / (estimateDfDv * 2**divisor); + bestValueImprovement = functionToMinimize(variableOneUpdate, variableTwoUpdate) ?? bestValue; + divisor += 1; + + if (divisor >= 32) { + deadEnd = true; + break; + } + } + + if (deadEnd) { + break; + } + + variableOne = variableOneUpdate; + variableTwo = variableTwoUpdate; + bestValue = bestValueImprovement; + } + + return [variableOne, variableTwo]; +} + export function calculateSimplePlaneChange(coordinates: OrbitalCoordinates, planet: Body, targetInclination: number, targetLongitudeOfAscendingNode: number, circularizeOrbit: boolean): SimplePlaneChange { const otherPlaneNormal = [ [Math.sin(targetLongitudeOfAscendingNode)*Math.sin(targetInclination)], @@ -807,8 +898,8 @@ export function findLambertSolutionsWithCorrectTime(lambertSolutions: LambertSol export function findCheapestIntercept(startingSituation: OrbitalCoordinates, targetSituation: OrbitalCoordinates, body: Body, extraTrueAnomaly: number, progressCallback?: ProgressCallbackFunction): Transfer | null { // Find acceptable intercept times // First number is intercept time, second number is true anomaly - let starts: [number, number][] = []; - let intercepts: [number, number][] = []; + let startTimes: number[] = []; + let endTimes: number[] = []; // Check if the starting orbit is stable0 let startingOrbitStable = false; @@ -828,95 +919,147 @@ export function findCheapestIntercept(startingSituation: OrbitalCoordinates, tar } } - let maxStartTime = 1e99; - let maxStartTrueAnomaly = 1e99; - let maxEndTime = 1e99; - let maxEndTrueAnomaly = 1e99; + let maxEndTime: number | null = null; + let maxStartTime: number | null = null; - if (interceptOrbitStable && startingOrbitStable) { - // If both orbits are stable, we'll check for three whole orbits of the largest orbit - let startingOrbitPeriod = getTimeBetweenTrueAnomalies(0, 2*Math.PI, startingSituation.orbit, body); - let targetOrbitPeriod = getTimeBetweenTrueAnomalies(0, 2*Math.PI, targetSituation.orbit, body); - - maxStartTime = 3 * Math.max(startingOrbitPeriod, targetOrbitPeriod); - maxEndTime = 4 * Math.max(startingOrbitPeriod, targetOrbitPeriod); - } else { - maxStartTime = 1e99; - maxStartTrueAnomaly = 1e99; - maxEndTime = 1e99; - maxEndTrueAnomaly = 1e99; - - if (!startingOrbitStable) { - maxStartTrueAnomaly = Math.abs(Math.acos((startingSituation.orbit.semiLatusRectum - body.sphereOfInfluence) / (body.sphereOfInfluence * startingSituation.orbit.eccentricity))); - maxStartTime = getTimeBetweenTrueAnomalies(startingSituation.trueAnomaly, maxStartTrueAnomaly, startingSituation.orbit, body); - } - - if (!interceptOrbitStable) { - maxEndTrueAnomaly = Math.abs(Math.acos((targetSituation.orbit.semiLatusRectum - body.sphereOfInfluence) / (body.sphereOfInfluence * targetSituation.orbit.eccentricity))); - maxEndTime = getTimeBetweenTrueAnomalies(targetSituation.trueAnomaly, maxEndTrueAnomaly, targetSituation.orbit, body); - } else { - // Set the max end time to the max start time plus two full orbits of the target - maxEndTime = maxStartTime + 4 * Math.PI * Math.sqrt((targetSituation.orbit.semiLatusRectum / (1 - targetSituation.orbit.eccentricity**2))**3 / body.gravitationalParameter); - } + if (!startingOrbitStable) { + let maxStartTrueAnomaly = Math.abs(Math.acos((startingSituation.orbit.semiLatusRectum - body.sphereOfInfluence) / (body.sphereOfInfluence * startingSituation.orbit.eccentricity))); + maxStartTime = getTimeBetweenTrueAnomalies(startingSituation.trueAnomaly, maxStartTrueAnomaly, startingSituation.orbit, body); + // If we're leaving the galaxy, set the max end time to the time it would take to orbit at the edge of the sphere of influence + maxEndTime = getOrbitalPeriod(body.sphereOfInfluence, body); } - let anomalyCounter = 1; - while (true) { - let possibleStart = startingSituation.trueAnomaly + anomalyCounter * 2 * Math.PI / 100; - if (possibleStart > maxStartTrueAnomaly) { - break; - } - let possibleStartTime = getTimeBetweenTrueAnomalies(startingSituation.trueAnomaly, possibleStart, startingSituation.orbit, body); - if (possibleStartTime > maxStartTime) { - break; - } - - starts.push([possibleStartTime, possibleStart]); - anomalyCounter += 1; + if (!interceptOrbitStable) { + let maxEndTrueAnomaly = Math.abs(Math.acos((targetSituation.orbit.semiLatusRectum - body.sphereOfInfluence) / (body.sphereOfInfluence * targetSituation.orbit.eccentricity))); + maxEndTime = Math.min(maxEndTime ?? 1e99, getTimeBetweenTrueAnomalies(targetSituation.trueAnomaly, maxEndTrueAnomaly, targetSituation.orbit, body)); } - anomalyCounter = 0; - while (true) { - let possibleEnd = targetSituation.trueAnomaly + anomalyCounter * 2 * Math.PI / 100.0; - if (possibleEnd > maxEndTrueAnomaly) { - break; - } - let possibleEndTime = getTimeBetweenTrueAnomalies(targetSituation.trueAnomaly, possibleEnd, targetSituation.orbit, body); - possibleEnd += extraTrueAnomaly; + if (startingOrbitStable && interceptOrbitStable) { + // If both orbits are stable, try for four of the longest orbit + let startingSemiMajor = startingSituation.orbit.semiLatusRectum / (1 - startingSituation.orbit.eccentricity**2); + let targetSemiMajor = targetSituation.orbit.semiLatusRectum / (1 - targetSituation.orbit.eccentricity**2); - if (possibleEndTime > maxEndTime) { - break; - } + let startingOrbitalPeriod = getOrbitalPeriod(startingSemiMajor, body); + let targetOrbitalPeriod = getOrbitalPeriod(targetSemiMajor, body); - intercepts.push([possibleEndTime, possibleEnd]); - anomalyCounter += 1; + maxStartTime = 4 * Math.max(startingOrbitalPeriod, targetOrbitalPeriod); + maxEndTime = 5 * Math.max(startingOrbitalPeriod, targetOrbitalPeriod); } - let pairs: [number, number, number, number][] = []; - starts.forEach(([startingTime, startingTrueAnomaly]) => { - intercepts.forEach(([endingTime, endingTrueAnomaly]) => { - if (endingTime > startingTime) { - pairs.push([startingTime, startingTrueAnomaly, endingTime, endingTrueAnomaly]); + // Max end time should have been set by now, but Typescript doesn't seem to know that + maxEndTime = maxEndTime ?? 0; + + // No point in having start times after the end times + if (!maxStartTime || maxStartTime > maxEndTime) { + maxStartTime = maxEndTime; + } + + // We'll try starting at 100 different start and end times + let startTimeStep = maxStartTime / 100; + let endTimeStep = maxEndTime / 100; + + for (var i = 0; i < 100; i++) { + startTimes.push(i * startTimeStep); + endTimes.push(i * endTimeStep); + } + + // Create pairs to test + let pairs: [number, number][] = []; + startTimes.forEach(startTime => { + endTimes.forEach(endTime => { + if (endTime > startTime) { + pairs.push([startTime, endTime]); } - }); + }) }); - let bestDeltaV: number | null = null; + const getTransfer = (startTime: number, endTime: number, backwards: boolean): Transfer | null => { + let startOrbitMeanAnomaly = startingSituation.meanAnomaly + startTime * startingSituation.meanAngularMotion; + let endOrbitMeanAnomaly = targetSituation.meanAnomaly + endTime * targetSituation.meanAngularMotion; + + let extraStartOrbits = Math.floor(startOrbitMeanAnomaly / (2 * Math.PI)); + let extraEndOrbits = Math.floor(endOrbitMeanAnomaly / (2 * Math.PI)); + + startOrbitMeanAnomaly = startOrbitMeanAnomaly % (2 * Math.PI); + endOrbitMeanAnomaly = endOrbitMeanAnomaly % (2 * Math.PI); + + let [_, startTrueAnomaly] = getEccentricAndTrueAnomalyFromMeanAnomaly(startOrbitMeanAnomaly, startingSituation.orbit.eccentricity); + let [__, endTrueAnomaly] = getEccentricAndTrueAnomalyFromMeanAnomaly(endOrbitMeanAnomaly, targetSituation.orbit.eccentricity); + + startTrueAnomaly += extraStartOrbits * 2 * Math.PI; + endTrueAnomaly += extraEndOrbits * 2 * Math.PI; + + let targetTransferTime = endTime - startTime; + let lambertSolutions = new LambertSolutions(startingSituation.orbit, startTrueAnomaly, targetSituation.orbit, endTrueAnomaly + extraTrueAnomaly, body, backwards); + + let transfer = findLambertSolutionsWithCorrectTime(lambertSolutions, targetTransferTime); + if (transfer) { + transfer.firstManoeuvre.time += startTime; + transfer.secondManoeuvre.time += startTime; + } + + return transfer; + } + + // Create the function we want to gradient descend + const getInterceptFunction = (backwards: boolean) => { + return function(startTime: number, endTime: number): number | null { + let transfer = getTransfer(startTime, endTime, backwards); + if (!transfer) { + return null; + } + + if (transfer.farthestPointDistance > body.sphereOfInfluence || transfer.closestPointDistance < body.closestSafeDistance) { + return null; + } + + return transfer.firstManoeuvre.totalDeltaV + transfer.secondManoeuvre.totalDeltaV; + } + }; + + // Start at each pair, and do a little gradient descending + let bestDeltaV: number | null = null;; let bestTransfer: Transfer | null = null; - pairs.forEach(([startingTime, startingTrueAnomaly, endingTime, endingTrueAnomaly], index) => { - let transferTime = endingTime - startingTime; - [true, false].forEach(goBackwards => { - let lambertSolutions = new LambertSolutions(startingSituation.orbit, startingTrueAnomaly, targetSituation.orbit, endingTrueAnomaly, body, goBackwards); - let transfer = findLambertSolutionsWithCorrectTime(lambertSolutions, transferTime); - if (transfer !== null && transfer.closestPointDistance > body.closestSafeDistance && transfer.farthestPointDistance < body.sphereOfInfluence) { - let totalDeltaV = transfer.firstManoeuvre.totalDeltaV + transfer.secondManoeuvre.totalDeltaV; - if (bestDeltaV == null || totalDeltaV < bestDeltaV) { - bestDeltaV = totalDeltaV; - bestTransfer = transfer; + pairs.forEach(([startTime, endTime], index) => { + // Try both forwards and backwards + [true, false].forEach(backwards => { + let interceptFunction = getInterceptFunction(backwards); - bestTransfer.firstManoeuvre.time += startingTime; - bestTransfer.secondManoeuvre.time += startingTime; + let foundStartTime = startTime; + let foundEndTime = endTime; + + let foundWorkingTransfer = false; + let trialCounter = 0; + while (!foundWorkingTransfer && trialCounter < 20) { + trialCounter++; + + let testTransfer = interceptFunction(foundStartTime, foundEndTime); + if (testTransfer) { + foundWorkingTransfer = true; + } else { + foundStartTime = Math.max(0,startTime + Math.random() * 0.5 * startTimeStep); + foundEndTime = Math.max(0,endTime + Math.random() * 0.5 * endTimeStep); + + if (foundEndTime < foundStartTime) { + foundEndTime += startTimeStep; + } + } + } + + if (!foundWorkingTransfer) { + return; + } + + let result = perform2dGradientDescent(interceptFunction, foundStartTime, foundEndTime, 0, 30); + if (result) { + let foundTransfer = getTransfer(result[0], result[1], backwards); + if (foundTransfer) { + let transferDeltaV = foundTransfer.firstManoeuvre.totalDeltaV + foundTransfer.secondManoeuvre.totalDeltaV; + if (bestDeltaV == null || transferDeltaV < bestDeltaV) { + bestDeltaV = transferDeltaV; + bestTransfer = foundTransfer; + } } } }); @@ -925,7 +1068,7 @@ export function findCheapestIntercept(startingSituation: OrbitalCoordinates, tar progressCallback(index+1, pairs.length, bestDeltaV, bestTransfer); } }); - + return bestTransfer; } @@ -1257,11 +1400,11 @@ export function findOrbitThroughInterpolation(ownCoordinates: OrbitalCoordinates } - const getPossibleTargetLocationFunction = (ownOrbit: Orbit, ownAltitude: number, targetAltitude: number, distanceToTarget: number): ((angle: number) => number[][]) => { + const getPossibleTargetVectors = (ownOrbit: Orbit, ownAltitude: number, targetAltitude: number, distanceToTarget: number): [number[][], number[][], number[][]] => { let ownTrueAnomaly = getTrueAnomalyFromAltitude(ownAltitude, ownOrbit.semiLatusRectum, ownOrbit.eccentricity, ownShipHeadedInwards); - let ownLocalX = interpolationParameters.firstOwnAltitude * Math.cos(ownTrueAnomaly); - let ownLocalY = interpolationParameters.firstOwnAltitude * Math.sin(ownTrueAnomaly); + let ownLocalX = ownAltitude * Math.cos(ownTrueAnomaly); + let ownLocalY = ownAltitude * Math.sin(ownTrueAnomaly); let ownDirection = normalizeVector( addVector( @@ -1270,8 +1413,8 @@ export function findOrbitThroughInterpolation(ownCoordinates: OrbitalCoordinates ) ); - let perpendicularOne = getRandomPerpendicularVector(ownDirection); - let perpendicularTwo = normalizeVector(vectorCrossProduct(ownDirection, perpendicularOne)); + let perpendicularInPlane = normalizeVector(vectorCrossProduct(ownDirection, ownOrbit.coordinateAxes[2])); + let perpendicularOutOfPlane = normalizeVector(vectorCrossProduct(ownDirection, perpendicularInPlane)); // Find phase angle to target let phaseAngle = Math.acos((ownAltitude**2 + targetAltitude**2 - distanceToTarget**2) / (2 * ownAltitude * targetAltitude)); @@ -1280,95 +1423,173 @@ export function findOrbitThroughInterpolation(ownCoordinates: OrbitalCoordinates let distanceAlongNormal = targetAltitude * Math.sin(phaseAngle); let root = multiplyMatrixWithScalar(distanceAlongDirection, ownDirection); + let vectorOne = multiplyMatrixWithScalar(distanceAlongNormal, perpendicularInPlane); + let vectorTwo = multiplyMatrixWithScalar(distanceAlongNormal, perpendicularOutOfPlane); - return function(angle: number) { - let addition = addVector( - multiplyMatrixWithScalar(distanceAlongNormal * Math.cos(angle), perpendicularOne), - multiplyMatrixWithScalar(distanceAlongNormal * Math.sin(angle), perpendicularTwo) - ); - - return addVector(root, addition); - } + return [root, vectorOne, vectorTwo]; } - let firstTargetFunction = getPossibleTargetLocationFunction(ownCoordinates.orbit, interpolationParameters.firstOwnAltitude, interpolationParameters.firstTargetAltitude, interpolationParameters.firstDistance); - let secondTargetFunction = getPossibleTargetLocationFunction(ownCoordinates.orbit, interpolationParameters.secondOwnAltitude, interpolationParameters.secondTargetAltitude, interpolationParameters.secondDistance); + const getFunctionAndDerivatives = (possiblePositionOne: [number[][], number[][], number[][]], possiblePositionTwo: [number[][], number[][], number[][]], altitudeOne: number, altitudeTwo: number, anomalyDifference: number): [(angleOne: number, angleTwo: number) => number, (angleOne: number, angleTwo: number) => number, (angleOne: number, angleTwo: number) => number] => { + let c1 = possiblePositionOne[0]; + let p11 = possiblePositionOne[1]; + let p12 = possiblePositionOne[2]; + + let c2 = possiblePositionTwo[0]; + let p21 = possiblePositionTwo[1]; + let p22 = possiblePositionTwo[2]; + + // Scale everything to avoid numerical explosions + let scaleFactor = 1 / getVectorMagnitude(c2); + + c1 = multiplyMatrixWithScalar(scaleFactor, c1); + p11 = multiplyMatrixWithScalar(scaleFactor, p11); + p12 = multiplyMatrixWithScalar(scaleFactor, p12); + c2 = multiplyMatrixWithScalar(scaleFactor, c2); + p21 = multiplyMatrixWithScalar(scaleFactor, p21); + p22 = multiplyMatrixWithScalar(scaleFactor, p22); + + let scaledAltitudeOne = altitudeOne * scaleFactor; + let scaledAltitudeTwo = altitudeTwo * scaleFactor; + + const constantPart = scaledAltitudeOne * scaledAltitudeTwo * Math.cos(anomalyDifference); + + const functionToMinimize = (angleOne: number, angleTwo: number) => { + let d1 = addVector(c1, addVector( + multiplyMatrixWithScalar(Math.cos(angleOne), p11), + multiplyMatrixWithScalar(Math.sin(angleOne), p12) + )); + + let d2 = addVector(c2, addVector( + multiplyMatrixWithScalar(Math.cos(angleTwo), p21), + multiplyMatrixWithScalar(Math.sin(angleTwo), p22) + )); + + return vectorDotProduct(d1, d2) - constantPart; + } + + const partialDerivativeAngleOne = (angleOne: number, angleTwo: number) => { + let d2 = addVector(c2, addVector( + multiplyMatrixWithScalar(Math.cos(angleTwo), p21), + multiplyMatrixWithScalar(Math.sin(angleTwo), p22) + )); + + let d1Dx = addVector( + multiplyMatrixWithScalar(-Math.sin(angleOne), p11), + multiplyMatrixWithScalar(Math.cos(angleOne), p12) + ); + + return vectorDotProduct(d1Dx, d2); + } + + const partialDerivativeAngleTwo = (angleOne: number, angleTwo: number) => { + let d1 = addVector(c1, addVector( + multiplyMatrixWithScalar(Math.cos(angleOne), p11), + multiplyMatrixWithScalar(Math.sin(angleOne), p12) + )); + + let d2Dy = addVector( + multiplyMatrixWithScalar(-Math.sin(angleTwo), p21), + multiplyMatrixWithScalar(Math.cos(angleTwo), p22) + ); + + return vectorDotProduct(d1, d2Dy); + } + + return [functionToMinimize, partialDerivativeAngleOne, partialDerivativeAngleTwo]; + } + + let firstPositionPossibleVectors = getPossibleTargetVectors(ownCoordinates.orbit, interpolationParameters.firstOwnAltitude, interpolationParameters.firstTargetAltitude, interpolationParameters.firstDistance); + let secondPositionPossibleVectors = getPossibleTargetVectors(ownCoordinates.orbit, interpolationParameters.secondOwnAltitude, interpolationParameters.secondTargetAltitude, interpolationParameters.secondDistance); // Since we know a lot about the target orbit, we can figure out how far it is supposed to be between the two points let firstTargetTrueAnomaly = getTrueAnomalyFromAltitude(interpolationParameters.firstTargetAltitude, targetSemiLatusRectum, targetEccentricity, targetHeadedInwards); let secondTargetTrueAnomaly = getTrueAnomalyFromAltitude(interpolationParameters.secondTargetAltitude, targetSemiLatusRectum, targetEccentricity, targetHeadedInwards); - // Use cosine rule - let targetMovement = Math.sqrt(interpolationParameters.firstTargetAltitude**2 + interpolationParameters.secondTargetAltitude**2 - 2 * interpolationParameters.firstTargetAltitude * interpolationParameters.secondTargetAltitude * Math.cos(secondTargetTrueAnomaly - firstTargetTrueAnomaly)); + let [minimizeFunction, partialDerivativeAngleOne, partialDerivativeAngleTwo] = getFunctionAndDerivatives(firstPositionPossibleVectors, secondPositionPossibleVectors, interpolationParameters.firstTargetAltitude, interpolationParameters.secondTargetAltitude, secondTargetTrueAnomaly - firstTargetTrueAnomaly); - // Now, we need to do multivariate optimization to find the two angles to give to the functions that make the distance between - // two target positions as equal to the target movement as possible - const functionToOptimize = (angleOne: number, angleTwo: number) => { - let firstTargetPosition = firstTargetFunction(angleOne); - let secondTargetPosition = secondTargetFunction(angleTwo); + const gradientDescent = (startingGuessAngleOne: number, startingGuessAngleTwo: number): [number, number, number] => { + let angleOne = startingGuessAngleOne; + let angleTwo = startingGuessAngleTwo; - return Math.abs(getVectorMagnitude(subtractVector(secondTargetPosition, firstTargetPosition)) - targetMovement); + let value = minimizeFunction(angleOne, angleTwo); + // Do up to 100 minimizing steps + for (var i = 0; i < 100; i++) { + let dfDx = partialDerivativeAngleOne(angleOne, angleTwo); + let dfDy = partialDerivativeAngleTwo(angleOne, angleTwo); + + let scaling = Math.sqrt(dfDx**2 + dfDy**2); + let directionX = -dfDy / scaling; + let directionY = -dfDx / scaling; + + let changeAlongDirection = directionX * dfDx + directionY * dfDy; + + let deadEnd = false; + let trialValue = value; + let trialAngleOne = angleOne; + let trialAngleTwo = angleTwo; + let divisor = 1; + + while (Math.abs(trialValue) >= Math.abs(value)) { + let stepAngleOne = - value * directionX / (changeAlongDirection * 2 ** divisor); + let stepAngleTwo = - value * directionY / (changeAlongDirection * 2 ** divisor); + trialAngleOne = angleOne + stepAngleOne; + trialAngleTwo = angleTwo + stepAngleTwo; + trialValue = minimizeFunction(trialAngleOne, trialAngleTwo); + divisor += 1; + + if (divisor >= 32) { + deadEnd = true; + break; + } + } + + if (deadEnd) { + return [angleOne, angleTwo, value]; + } + + value = trialValue; + angleOne = trialAngleOne; + angleTwo = trialAngleTwo; + } + + return [angleOne, angleTwo, value]; } let bestAngleOne = 0; let bestAngleTwo = 0; - let bestValue = 1e99; + let bestValue = null; - for (var i = 0; i < 100; i++) { - for (var j = 0; j < 100; j++) { - let angleOne = 2 * Math.PI * i / 100; - let angleTwo = 2 * Math.PI * j / 100; + // Try to start at a bunch of different positions and do gradient descent from there + for (var i = 0; i < 500; i++) { + for (var j = 0; j < 500; j++) { + let startingAngleOne = i * 2 * Math.PI / 500.0; + let startingAngleTwo = j * 2 * Math.PI / 500.0 + let startingValue = minimizeFunction(startingAngleOne, startingAngleTwo); - let value = functionToOptimize(angleOne, angleTwo); - if (value < bestValue) { - bestAngleOne = angleOne; - bestAngleTwo = angleTwo; - bestValue = value; + if (bestValue == null || Math.abs(startingValue) < Math.abs(bestValue)) { + bestAngleOne = startingAngleOne; + bestAngleTwo = startingAngleTwo; + bestValue = startingValue; } } } - // Do some gradient descent on the best values found so far - while (bestValue > 0.5) { - let estimateDfDx = (functionToOptimize(bestAngleOne+0.0000001, bestAngleTwo) - bestValue) / 0.0000001; - let estimateDfDy = (functionToOptimize(bestAngleOne, bestAngleTwo+0.0000001) - bestValue) / 0.0000001; + [bestAngleOne, bestAngleOne, bestValue] = gradientDescent(bestAngleOne, bestAngleTwo); - let slopeNormal = vectorCrossProduct( - [[1], [0], [estimateDfDx]], - [[0], [1], [estimateDfDy]] - ); - - let direction = normalizeVector([[slopeNormal[0][0]], [slopeNormal[1][0]], [0]]); - let estimateDfDv = (functionToOptimize(bestAngleOne + direction[0][0]*0.0001, bestAngleTwo + direction[1][0]*0.0001) - bestValue) / 0.0001; - - let bestValueImprovement = bestValue; - let bestAngleOneUpdate = bestAngleOne; - let bestAngleTwoUpdate = bestAngleTwo; - let divisor = 0; - let deadEnd = false; - - while (bestValueImprovement >= bestValue) { - bestAngleOneUpdate = bestAngleOne - bestValue * direction[0][0] / (estimateDfDv * 2**divisor); - bestAngleTwoUpdate = bestAngleTwo - bestValue * direction[1][0] / (estimateDfDv * 2**divisor); - bestValueImprovement = functionToOptimize(bestAngleOneUpdate, bestAngleTwoUpdate); - divisor += 1; - - if (divisor >= 32) { - deadEnd = true; - break; - } - } - - if (deadEnd) { - break; - } - bestAngleOne = bestAngleOneUpdate; - bestAngleTwo = bestAngleTwoUpdate; - bestValue = bestValueImprovement; - } - - let bestFirstPosition = firstTargetFunction(bestAngleOne); - let bestSecondPosition = secondTargetFunction(bestAngleTwo); + let bestFirstPosition = addVector( + firstPositionPossibleVectors[0], + addVector( + multiplyMatrixWithScalar(Math.cos(bestAngleOne), firstPositionPossibleVectors[1]), + multiplyMatrixWithScalar(Math.sin(bestAngleOne), firstPositionPossibleVectors[2]) + ) + ); + let bestSecondPosition = addVector( + secondPositionPossibleVectors[0], + addVector( + multiplyMatrixWithScalar(Math.cos(bestAngleTwo), secondPositionPossibleVectors[1]), + multiplyMatrixWithScalar(Math.sin(bestAngleTwo), secondPositionPossibleVectors[2]) + ) + ); if (bestFirstPosition == null || bestSecondPosition == null) { return null; diff --git a/src/gui/worker.ts b/src/gui/worker.ts index 8dadc8b..aad183e 100644 --- a/src/gui/worker.ts +++ b/src/gui/worker.ts @@ -46,7 +46,7 @@ export interface LandingProgressMessage { ctx.addEventListener("message", event => { const progressCallback = (numberChecked: number, totalNumber: number, bestDeltaV: number | null, bestTransfer: Transfer | null) => { - if (numberChecked % 100 == 0) { + if (numberChecked % 1 == 0) { let percentDone = numberChecked * 100 / totalNumber; let message: ProgressMessage = { type: "ProgressMessage",