Алгоритм планирования встречи между двумя космическими кораблями

Question

Я пытаюсь выяснить алгоритм установки рандеву между двумя космическими кораблями.

Нет гравитации или сопротивления. Оба космических корабля имеют позицию и скорость на старте. Космический корабль B продолжает свой курс без ускорения, поэтому космическому кораблю A необходимо ускориться, чтобы сократить расстояние между ними и затем согласовать скорости, когда он достигнет положения космического корабля B.

Космический корабль может мгновенно изменить свое направление тяга, но может использовать только максимальное ускорение или вообще не использовать ускорение. Я также хочу ограничить разность скоростей между космическими кораблями во время маневра.

Мне бы хотелось, чтобы выходные данные были в виде ряда ветвей траектории, то есть: leg1: направление ускорения x в течение t1 секунд, leg2: уклон в течение t2 секунд, leg3: ускорение в направлении y в течение t3 секунд.

Мне не нужно оптимальное решение, но я бы хотел, чтобы оно «выглядело правильно».

Я пытался создать импульс, чтобы выровнять скорости и добавить его к импульсу для движения к космическому кораблю B, но даже если космический корабль A достигает правильной скорости, он не может достичь позиции цели. Я попробовал импульсы самостоятельно, и они, кажется, работают должным образом, поэтому я предполагаю, что проблема в том, как я их складываю. Я не знаю, правильно ли я это реализую или этот подход просто не сработает. Я надеюсь, что кто-то с более сильными навыками математики и физики может просветить меня.

Вот код, который я использую:

// var velocityAdjustmentTime =  (int)Math.Sqrt(2 * velocityDelta.Length / tp.Acceleration);
            var velocityAdjustmentTime = (int)(velocityDelta.Length / tp.Acceleration);

            var velocityAdjustVector = velocityDelta;
            velocityAdjustVector.Normalize();
            velocityAdjustVector *= tp.Acceleration;

            var targetAccelerationDisplacement = new Vector3D(0, 0, 0); // TODO: Replace this with proper values.

            Vector3D newPosition;
            Vector3D newVelocity;
            Vector3D targetNewPosition;

            // Check if the formation and the target already have a paralell course with the same velocity.
            if (velocityAdjustmentTime > 0)
            {
                // If not, calculate the position and velocity after the velocity has been aligned.
                newPosition = tp.StartPosition + (tp.StartVelocity * velocityAdjustmentTime) + ((velocityAdjustVector * velocityAdjustmentTime * velocityAdjustmentTime) / 2);
                newVelocity = tp.StartVelocity + velocityAdjustVector * velocityAdjustmentTime;
                targetNewPosition = tp.TargetStartPosition + (tp.TargetStartVelocity * velocityAdjustmentTime) + targetAccelerationDisplacement;
            }

            else
            {
                // Else, new and old is the same.
                newPosition = tp.StartPosition;
                newVelocity = tp.StartVelocity;
                targetNewPosition = tp.TargetStartPosition;
            }

            // Get the new direction from the position after velocity change.
            var newDirection = targetNewPosition - newPosition;


            // Changing this value moves the end position closer to the target. Thought it would be newdirection length, but then it doesn't reach the target.
            var length = newDirection.Length;

            // I don't think this value matters.
            var speed = (int)(cruiseSpeed);

            var legTimes = CalculateAccIdleDecLegs(tp.Acceleration, length, speed);

            // Sets how much of the velocity change happens on the acceleration or deceleration legs.
            var velFactorAcc = 1;
            var velFactorDec = 1 - velFactorAcc;

            // Make the acceleration vector.
            accelerationVector = newDirection;
            accelerationVector.Normalize();
            accelerationVector *= legTimes[0] * tp.Acceleration;

            accelerationVector += velocityDelta * velFactorAcc;



            accelerationTime = (int)(accelerationVector.Length / tp.Acceleration);

            accelerationVector.Normalize();
            accelerationVector *= tp.Acceleration;


            // Make the acceleration order.
            accelerationLeg.Acceleration = accelerationVector;
            accelerationLeg.Duration = accelerationTime;

            // Make the deceleration vector.
            decelerationVector = newDirection;
            decelerationVector.Negate();
            decelerationVector.Normalize();
            decelerationVector *= legTimes[2] * tp.Acceleration;

            decelerationVector += velocityDelta * velFactorDec;

            decelerationTime = (int)(decelerationVector.Length / tp.Acceleration);

            decelerationVector.Normalize();
            decelerationVector *= tp.Acceleration;

            // And deceleration order.
            decelerationLeg.Acceleration = decelerationVector;
            decelerationLeg.Duration = decelerationTime;

            // Add the orders to the list.
            trajectory.Add(accelerationLeg);

            // Check if there is an idle leg in the middle...
            if (legTimes[1] > 0)
            {
                // ... if so, make the order and add it to the list.
                idleLeg.Duration = legTimes[1];

                trajectory.Add(idleLeg);
            }

            // Add the deceleration order.
            trajectory.Add(decelerationLeg);

И функция для вычисления подходов:

private static int[] CalculateAccIdleDecLegs(double acceleration, double distance, int cruiseSpeed)
    {
        int[] legDurations = new int[3];
        int accelerationTime;
        int idleTime;
        int decelerationTime;

        // Calculate the max speed it's possible to accelerate before deceleration needs to begin.
        var topSpeed = Math.Sqrt(acceleration * distance);

        // If the cruise speed is higher than or equal to the possible top speed, the formation should accelerate to top speed and then decelerate.
        if (cruiseSpeed >= topSpeed)
        {
            // Get the time to accelerate to the max velocity.
            accelerationTime = (int)((topSpeed) / acceleration);

            // Idle time is zero.
            idleTime = 0;

            // Get the deceleration time.
            decelerationTime = (int)(topSpeed / acceleration);
        }

        // Else, the formation should accelerate to max velocity and then coast until it starts decelerating.
        else
        {
            // Find the acceleration time.
            accelerationTime = (int)((cruiseSpeed) / acceleration);

            // Get the deceleration time.
            decelerationTime = (int)(cruiseSpeed / acceleration);

            // Calculate the distance traveled while accelerating.
            var accelerationDistance = 0.5 * acceleration * accelerationTime * accelerationTime;

            // Calculate the distance traveled while decelerating.
            var decelerationDistance = 0.5 * acceleration * decelerationTime * decelerationTime;

            // Add them together.
            var thrustDistance = accelerationDistance + decelerationDistance;

            // Find the idle distance.
            var idleDistance = distance - thrustDistance;

            // And the time to idle.
            idleTime = (int)(idleDistance / cruiseSpeed);
        }

        legDurations[0] = accelerationTime;
        legDurations[1] = idleTime;
        legDurations[2] = decelerationTime;

        return legDurations;
    }

Answer 1

НОВАЯ ВЕРСИЯ:

Предположим, у вас есть начальные позиции и скорости xA0, vA0 и xB0, vB0 космических кораблей A и B соответственно. Как вы сказали, B движется без ускорения и с постоянной скоростью vB0. Поэтому он равномерно движется по прямой. Его движение описывается как: xB = xB0 + t*vB0. Космический корабль А может включать и выключать ускорение постоянной величины a0, но может менять направление по своему усмотрению. Скорость A не должна превышать определенного значения v_max > 0.

Поскольку космический корабль B движется равномерно по прямой с постоянной скоростью vB0, он фактически определяет инерциальную систему координат. Другими словами, если мы переведем исходную систему координат и прикрепим ее к B, новая система движется с постоянной скоростью вдоль прямой линии и поэтому также является инерционной. Преобразование галилеево, поэтому можно определить следующее изменение координат (в обоих направлениях)

y = x - xB0 - t*vB0
u = v - vB0

x = y + xB0 + t*vB0
v = u + vB0

, в частности, для B для любого момента времени t мы получим

yB = xB - xB0 - t*vB0 = xB0 + t*vB0 - xB0 - t*vB0  = 0``

В момент t=0,

yA0 = xA0 - xB0  

uA0 = vA0 - vB0

Мы стремимся спроектировать элемент управления в этой новой системе координат, и они переместят его обратно в исходную. Итак, давайте переключим координаты:

y = x - xB
u = v - vB0

Итак, в этой новой инерциальной системе координат мы решаем задачу теории управления, и для создания хорошего управления мы будем использовать в качестве функции Ляпунова (функцию, которая позволяет нам чтобы гарантировать определенное устойчивое поведение и разработать правильное выражение для ускорения a), величина скорости в квадрате L = norm(u)^2. Мы хотим спроектировать ускорение a, чтобы функция Ляпунова в начальной фазе движения монотонно и неуклонно уменьшалась, а скорость переориентировалась соответствующим образом.

Определить единичный вектор

L_unit = cross(x0A - xB0, v0A - vB0) / norm(cross(x0A - xB0, v0A - vB0))

Пусть в системе координат, связанной с B, движение A удовлетворяет системе обыкновенных дифференциальных уравнений (эти уравнения в обеих системах являются уравнениями Ньютона, поскольку обе системы являются инерционными):

dy/dt = u

du/dt = - a0 * (u - cross(L_unit, u)) / norm(u - cross(L_unit, u))

Другими словами, ускорение установлено на

a = - a0 * (u - cross(L_unit, u)) / norm(u - cross(L_unit, u))

Обратите внимание, что norm(a) = a0. Поскольку векторы u и cross(L_unit, u) являются ортогональными и равными по величине (просто cross(L_unit, u) - это поворот вектора на * девяносто градусов *), знаменатель упрощается до

norm(u - cross(L_unit, u)) = sqrt( norm(u - cross(L_unit, u))^2 ) 
                           = sqrt( norm(u)^2 + norm(L_unit, u)^2 )
                           = sqrt( norm(u)^2 + norm(L_unit)^2*norm(u)^2 )
                           = sqrt( norm(u)^2 + norm(u)^2)
                           = sqrt(2) * norm(u)

Таким образом, система дифференциальные уравнения упрощаются до

dy/dt = u

du/dt = -(a0/sqrt(2)) * u/norm(u) + (a0/sqrt(2)) * cross(L_unit, u)) / norm(u)

Система спроектирована таким образом, что A всегда перемещается в плоскости, проходящей через начало координат B и перпендикулярно вектору L_unit.

Поскольку u и cross(L_unit, u) перпендикулярны, их точечное произведение равно 0, что позволяет нам вычислять производную по времени функции Ляпунова вдоль решений для системы выше (u' означает транспонирование column-vector u):

d/dt( L ) = d/dt( norm(u)^2 ) = d/dt( u' * u ) = u' * du/dt 
          = u' * (  -(a0/sqrt(2)) * u/norm(u) 
            + (a0/sqrt(2)) * cross(L_unit, u)) / norm(u) )
          = -(a0/sqrt(2)) * u'*u / norm(u) 
            + (a0/sqrt(2)) * u'*cross(L_unit, u)) / norm(u) 
          = -(a0/sqrt(2)) * norm(u)^2 / norm(u)
          = -(a0/sqrt(2)) * norm(u)
          = - (a0/sqrt(2)) * sqrt(L)

d/dt( L ) = -(a0/sqrt(2)) * sqrt(L) < 0

, что означает, что norm(u) уменьшается с течением времени до 0, как требуется.

Система дифференциальных уравнений, которая управляет движением, нелинейный и явно не разрешимый, поэтому мы должны интегрировать его численно. Чтобы сделать это, я выбрал геометрический метод интегратора c, в котором система разбита на две разрешимые для ясности системы, решения которых объединены вместе, чтобы дать (очень хорошее приближение) решение исходной системы. Это следующие системы:

dy/dt = u / 2

du/dt = -(a0/sqrt(2)) u / norm(u)

и

dy/dt = u / 2

du/dt = (a0/sqrt(2)) cross(L_unit, u) / norm(u)

Первоначально вторая система является нелинейной, однако после вычисления:

d/dt(norm(u)*2) = d/dt (dot(u, u)) = 2 * dot(u, du/dt) 
                = 2 * dot(u, (a0/sqrt(2)) * cross(L_unit , u))
                = 2 * (a0/sqrt(2)) * dot(u, cross(L_unit , u)) 
                = 0

мы заключаем, что во время При движении, определяемом этой системой, величина скорости постоянна, т.е. norm(u) = norm(u0), где u0 = u(0). Таким образом, системы, вместе со своими решениями, теперь выглядят так:

First system:

dy/dt = u / 2

du/dt = -(a0/sqrt(2)) u / norm(u)

Solution:
y(t) = y0  +  h * u0/2  -  t^2 * a0 * u0 / (4*sqrt(2)*norm(u0));
u(t) = u - t * a0 * u0 / (sqrt(2)*norm(u0));

и

Second system:

dy/dt = u / 2

du/dt = (a0/sqrt(2)// norm(u0)) cross(L_unit, u) 

Solution:
y(t) = y0 + (sqrt(2)*norm(u0)/a0) *( cross(L_unit, u0)
          + sin( t * a0/(sqrt(2)*norm(u0)) ) * u0  
          - cos( t  *a0/(sqrt(2)*norm(u0)) ) * cross(L_unit, u0) )     
u(t) = cos( t  *a0/(sqrt(2)*norm(u0)) ) * u0  
       + sin( t  *a0/(sqrt(2)*norm(u0)) ) * cross(L_unit, u0)

Решение исходной системы можно аппроксимировать следующим образом. Выберите шаг по времени h. Затем, если в момент времени t положение и скорость космического корабля были вычислены как y, u, положение и скорость обновленного космического корабля в момент времени t + h можно рассчитать, сначала позволив кораблю двигаться по решению второй системы, начиная с y, u для времени h/2, затем двигайтесь вдоль решения первой системы для времени h и затем двигайтесь вдоль решения второй системы для времени h/2.

function main()
    h = 0.3;
    a0 = 0.1;
    u_max = .8; % velocity threshold
    xB0 = [0; 0; 0];
    vB0 = [-1; 2; 0];
    xA0 = [ 7; 12; 0] + xB0;
    vA0 = [1; 5; 0]/7;
    %vA0 =  [2; -1; 0];
    L_unit = cross(xA0 - xB0, vA0 - vB0);
    L_unit =  L_unit / norm(L_unit);
    t = 0;
    xB = xB0;
    x = xA0;
    v = vA0;
    hold on
    grid on
    %axis([-200 20 -100 350])
    plot(0, 0, 'bo')

    % STEP 0 (the motion as described in the text above):
    n = floor(sqrt(2)*norm(vA0 - vB0)/(a0*h));
    for i=1:n
        [t, x, v, xB] = E(t, x, v, xB, vB0, a0, L_unit, h);
        plot(x(1), x(2), 'ro');
        plot(xB(1), xB(2), 'bo');
        pause(0.1)
    end
    u = v - vB0;
    norm_u = norm(u);
    % short additional decceleration so that A attains velocity v = vB0 
    t0 = t + norm_u/a0;    
    n = floor((t0 - t)/h); 
    a = - a0 * u / norm_u;    
    for i=1:n
        [t, x, v, xB] = ET(t, x, v, xB, vB0, a, h);
        plot(x(1), x(2), 'ro');
        plot(xB(1), xB(2), 'bo');
        pause(0.1)
    end
    [t, x, v, xB] = ET(t, x, v, xB, vB0, a, t0-t);
    plot(x(1), x(2), 'ro');
    plot(xB(1), xB(2), 'bo');
    pause(0.1)

    % STEP 1 (uniform acceleration of magnitude a0):
    v = vB0; 
    a = x-xB;
    norm_y0 = norm(a);
    a = - a0 * a / norm_y0;
    %t2 = t1 + sqrt( norm_y/a0 ); 
    accel_time = min( u_max/a0, sqrt( norm_y0/a0 ) );
    t1 = t0 + accel_time;
    n = floor((t1 - t0)/h);     
    for i=1:n
       [t, x, v, xB] = ET(t, x, v, xB, vB0, a, h);
       plot(x(1), x(2), 'bo');
       plot(xB(1), xB(2), 'ro');
       pause(0.1)
    end 
    [t, x, v, xB] = ET(t, x, v, xB, vB0, a, t1-t);
    plot(x(1), x(2), 'bo');
    plot(xB(1), xB(2), 'ro');
    pause(0.1)

    % STEP 2 (uniform straight-line motion): 
    norm_y1 = norm(x-xB);
    norm_y12 = max(0, norm_y0 - 2*(norm_y0 - norm_y1));
    t12 = norm_y12 / norm(v-vB0)
    t = t + t12
    n12 = floor(t12/h)
    for i=1:n12
        x = x + h*v; 
        xB = xB + h*vB0;
        plot(x(1), x(2), 'ro');
        plot(xB(1), xB(2), 'bo');
        pause(0.1)
    end
    x = x + (t12-n12*h)*v; 
    xB = xB + (t12-n12*h)*vB0;
    plot(x(1), x(2), 'ro');
    plot(xB(1), xB(2), 'bo');
    pause(0.1)

    % STEP 3 (uniform deceleration of magnitude a0, symmetric to STEP 1): 
    a = -a;
    for i=1:n % t2 + (t2-t1)
       [t, x, v, xB] = ET(t, x, v, xB, vB0, a, h);
       plot(x(1), x(2), 'bo');
       plot(xB(1), xB(2), 'ro');
       pause(0.1)
    end
    [t, x, v, xB] = ET(t, x, v, xB, vB0, a, t0+t12+2*accel_time-t);
    plot(x(1), x(2), 'bo');
    plot(xB(1), xB(2), 'ro');
    pause(0.1)
    norm(x-xB)
    norm(v-vB0)
end

Вот дополнительные функции, которые используются в главном коде выше:

% change of coordinates from world coordinates x, v 
% to coordinates y, u from spaship B's point of view:
function [y, u] = change(x, v, xB, vB0)
    y = x - xB;
    u = v - vB0;
end

% inverse chage of coordinates from y, u to x, v
function [x, v] = inv_change(y, u, xB, vB0)
    x = y + xB;
    v = u + vB0;
end

% solution to the second system of differential equations for a step h:
function [y_out, u_out] = R(y, u, a0, L_unit, h)
   omega = a0 / (sqrt(2) * norm(u));
   L_x_u = cross(L_unit, u);
   cos_omega_h = cos(omega*h);
   sin_omega_h = sin(omega*h);
   omega = 2*omega;
   y_out = y + (L_x_u ...
       + sin_omega_h * u  -  cos_omega_h * L_x_u) / omega;      
   u_out = cos_omega_h * u  +  sin_omega_h * L_x_u;
end

% solution to the first system of differential equations for a step h:
function [y_out, u_out] = T(y, u, a0, h)
    sqrt_2 = sqrt(2);
    u_unit = u / norm(u);  
    y_out = y  +  h * u/2  -  h^2 * a0 * u_unit/ (4*sqrt_2);
    u_out = u - h * a0 * u_unit / sqrt_2;
end

% approximate solution of the original system of differential equations for step h
% i.e. the sum of furst and second systems of differential equations:
function [t_out, x_out, v_out, xB_out] = E(t, x, v, xB, vB0, a0, L_unit, h)
   t_out = t + h;
   [y, u] = change(x, v, xB, vB0);
   [y, u] = R(y, u, a0, L_unit, h/2);
   [y, u] = T(y, u, a0, h);
   [y, u] = R(y, u, a0, L_unit, h/2);
   xB_out = xB + h*vB0;
   [x_out, v_out] = inv_change(y, u, xB_out, vB0);  
end

% straight-line motion with constant acceleration: 
function [t_out, x_out, v_out, xB_out] = ET(t, x, v, xB, vB0, a, h)
    t_out = t + h;
    [y, u] = change(x, v, xB, vB0);
    y = y  +  h * u  +  h^2 * a / 2;
    u = u + h * a;
    xB_out = xB + h*vB0;
    [x_out, v_out] = inv_change(y, u, xB_out, vB0);
end

СТАРШАЯ ВЕРСИЯ:

Я разработал два моделей. Обе модели изначально описываются в инерциальной системе отсчёта B y, u (см. Мои предыдущие ответы), а затем координаты преобразуются в исходные x, v. Я разработал управление на основе функции norm(u)^2 как функцию Ляпунова, так что на первом шаге алгоритма ускорение спроектировано так, чтобы функция Ляпунова norm(u)^2 неуклонно уменьшалась. В первой версии скорость уменьшения равна квадратичной c, но модель легче интегрировать, в то время как во второй версии скорость снижения является экспоненциальной, но модель требует интеграции Рунге-Кутты. И я не совсем настроил это. Я думаю, что Версия 1 должна выглядеть хорошо.

Взять L_unit = cross(y0, u0) / norm(cross(y0, u0)).

Версия 1: Модель:

dy/dt = y
du/dt = - a0 * (u + cross(L_unit, u)) / norm(u + cross(L_unit, u))
      = - a0 * (u + cross(L_unit, u)) / (norm(u)*sqrt(1 + norm(L_unit)^2))
      = - a0 * (u + cross(L_unit, u)) / (sqrt(2) * norm(u))

Чтобы интегрировать ее, разбейте ее на пару систем:

dy/dt = y
du/dt = - a0 * u / norm(u)

dy/dt = y
du/dt = - a0 * cross(L_unit, u) / norm(u0)  (see previous answers)

и интегрируйте их один за другим для небольших приращений h временных интервалов, а затем последовательно go назад и вперед между этими двумя системами. Я экспериментировал с некоторым кодом Matlab:

function main()
    h = 0.3;
    a0 = 0.1;
    xB0 = [0; 0; 0];
    vB0 = [-1; 2; 0];
    xA0 = [ 7; 12; 0] + xB0;
    vA0 =  [2; -1; 0];
    L_unit = cross(xA0 - xB0, vA0 - vB0);
    L_unit =  L_unit / norm(L_unit);
    t = 0;
    xB = xB0;
    x = xA0;
    v = vA0;
    hold on
    grid on
    %axis([-200 20 -100 350])
    plot(0, 0, 'bo')
    n = floor(2*norm(v - vB0)/(h*a0));
    for i=1:n
        [t, x, v, xB] = R(t, x, v, xB, vB0, a0, L_unit, h/2);
        a = - a0 * (v - vB0) / norm(v - vB0);
        [t, x, v, xB] = T(t, x, v, xB, vB0, a, h/2);
        plot(x(1), x(2), 'ro');
        plot(xB(1), xB(2), 'bo');
        pause(0.1)
    end
    t1 = t + norm(v - vB0)/a0;
    n = floor((t1 - t)/h); 
    a = - a0 * (v - vB0) / norm(v - vB0);
    for i=1:n
        [t, x, v, xB] = T(t, x, v, xB, vB0, a, h);
        plot(x(1), x(2), 'ro');
        plot(xB(1), xB(2), 'bo');
        pause(0.1)
    end
    [t, x, v, xB] = T(t, x, v, xB, vB0, a, t1-t);
    plot(x(1), x(2), 'ro');
    plot(xB(1), xB(2), 'bo');
    pause(0.1)
    t2 = t1 + sqrt( norm(x - xB)/a0 ); 
    n = floor((t2 - t1)/h);
    a = - a0 * (x - xB) / norm(x - xB);
    v = vB0;
    for i=1:n
        [t, x, v, xB] = T(t, x, v, xB, vB0, a, h);
        plot(x(1), x(2), 'ro');
        plot(xB(1), xB(2), 'bo');
        pause(0.1) 
    end
    [t, x, v, xB] = T(t, x, v, xB, vB0, a, t2-t);
    plot(x(1), x(2), 'ro');
    plot(xB(1), xB(2), 'bo');
    pause(0.1)
    for i=1:n % t2 + (t2-t1)
       [t, x, v, xB] = T(t, x, v, xB, vB0, -a, h);
       plot(x(1), x(2), 'ro');
       plot(xB(1), xB(2), 'bo');
       pause(0.1) 
    end
    [t, x, v, xB] = T(t, x, v, xB, vB0, -a, 2*t2 - t1 -t);
    plot(x(1), x(2), 'ro');
    plot(xB(1), xB(2), 'bo');
    pause(0.1)
end

, где соответствующие функции:

function [t_out, y_out, u_out] = R1(t, y, u, a0, L_unit, h)
   t_out = t + h;
   norm_u = norm(u);
   R = norm_u^2 / a0;
   cos_omega_h = cos(a0 * h / norm_u);
   sin_omega_h = sin(a0 * h / norm_u);
   u_unit = u / norm_u;
   y_out = y + R * cross(L_unit, u_unit) ...
           + R * sin_omega_h * u_unit ...
           - R * cos_omega_h * cross(L_unit, u_unit);
   u_out = norm_u * sin_omega_h * cross(L_unit, u_unit) ...
           + norm_u * cos_omega_h * u_unit;
end

function [t_out, x_out, v_out, xB_out] = R(t, x, v, xB, vB0, a0, L_unit, h)
    [t_out, y_out, u_out] = R1(t, x - xB, v - vB0, a0, L_unit, h);
    xB_out = xB + h * vB0;
    x_out = y_out + xB_out;
    v_out = u_out + vB0;
end

function [t_out, y_out, u_out] = T1(t, y, u, a, h)
    t_out = t + h;
    u_out = u + h * a;
    y_out = y + h * u + h^2 * a / 2; 
end

function [t_out, x_out, v_out, xB_out] = T(t, x, v, xB, vB0, a, h)
    [t_out, y_out, u_out] = T1(t, x - xB, v - vB0, a, h);
    xB_out = xB + h * vB0;
    x_out = y_out + xB_out;
    v_out = u_out + vB0;
end

Версия 2: Модель:

0 < k0 < 2 * a0 / norm(u0)  

dy/dt = y
du/dt = - k0 * u / 2 + sqrt(a0^2 - k0^2 * norm_u^2 / 4) * cross(L_unit, u/norm_u);

Код Matlab:

function main()
    h = 0.3;
    a0 = 0.1;
    xB0 = [0; 0; 0];
    vB0 = [-1; 2; 0];
    xA0 = [ 7; 12; 0] + xB0;
    vA0 =  [2; -1; 0];
    k0 = a0/norm(vA0-vB0);
    L_unit = cross(xA0 - xB0, vA0 - vB0);
    L_unit =  L_unit / norm(L_unit);
    t = 0;
    xB = xB0;
    x = xA0;
    v = vA0;
    hold on
    grid on
    %axis([-200 20 -100 350])
    plot(0, 0, 'bo')
    n = floor(2*norm(v - vB0)/(h*a0)); % this needs to be improved 
    for i=1:n
        [t, x, v, xB] = F_step(t, x, v, xB, vB0, a0, L_unit, k0, h);
        plot(x(1), x(2), 'ro');
        plot(xB(1), xB(2), 'bo');
        pause(0.1)
    end
    t1 = t + norm(v - vB0)/a0;
    n = floor((t1 - t)/h); 
    a = - a0 * (v - vB0) / norm(v - vB0);
    for i=1:n
        [t, x, v, xB] = T(t, x, v, xB, vB0, a, h);
        plot(x(1), x(2), 'ro');
        plot(xB(1), xB(2), 'bo');
        pause(0.1)
    end
    [t, x, v, xB] = T(t, x, v, xB, vB0, a, t1-t);
    plot(x(1), x(2), 'ro');
    plot(xB(1), xB(2), 'bo');
    pause(0.1)
    t2 = t1 + sqrt( norm(x - xB)/a0 ); 
    n = floor((t2 - t1)/h);
    a = - a0 * (x - xB) / norm(x - xB);
    v = vB0;
    for i=1:n
        [t, x, v, xB] = T(t, x, v, xB, vB0, a, h);
        plot(x(1), x(2), 'ro');
        plot(xB(1), xB(2), 'bo');
        pause(0.1) 
    end
    [t, x, v, xB] = T(t, x, v, xB, vB0, a, t2-t);
    plot(x(1), x(2), 'ro');
    plot(xB(1), xB(2), 'bo');
    pause(0.1)
    for i=1:n % t2 + (t2-t1)
       [t, x, v, xB] = T(t, x, v, xB, vB0, -a, h);
       plot(x(1), x(2), 'ro');
       plot(xB(1), xB(2), 'bo');
       pause(0.1) 
    end
    [t, x, v, xB] = T(t, x, v, xB, vB0, -a, 2*t2 - t1 -t);
    plot(x(1), x(2), 'ro');
    plot(xB(1), xB(2), 'bo');
    pause(0.1)
end

, где соответствующие функции:

function [dydt, dudt] = F1(u, a0, L_unit, k0)
    norm_u = norm(u);
    dydt = u;
    dudt = - k0 * u / 2 + sqrt(a0^2 - k0^2 * norm_u^2/4) * cross(L_unit, u/norm_u);
end

function [t_out, y_out, u_out] = F1_step(t, y, u, a0, L_unit, k0, h)
    t_out = t + h;
    [z1, w1] = F1(u, a0, L_unit, k0);
    [z2, w2] = F1(u + h * w1/2, a0, L_unit, k0);
    [z3, w3] = F1(u + h * w2/2, a0, L_unit, k0);
    [z4, w4] = F1(u + h * w3, a0, L_unit, k0);
    y_out = y + h*(z1 + 2*z2 + 2*z3 + z4)/6;
    u_out = u + h*(w1 + 2*w2 + 2*w3 + w4)/6;
end

function [t_out, x_out, v_out, xB_out] = F_step(t, x, v, xB, vB0, a0, L_unit, k0, h)
    [t_out, x_out, v_out] = F1_step(t, x-xB, v-vB0, a0, L_unit, k0, h);
    xB_out = xB + h * vB0;
    x_out = x_out + xB_out;
    v_out = v_out + vB0;
end

function [t_out, y_out, u_out] = T1(t, y, u, a, h)
    t_out = t + h;
    u_out = u + h * a;
    y_out = y + h * u + h^2 * a / 2; 
end

function [t_out, x_out, v_out, xB_out] = T(t, x, v, xB, vB0, a, h)
    [t_out, y_out, u_out] = T1(t, x - xB, v - vB0, a, h);
    xB_out = xB + h * vB0;
    x_out = y_out + xB_out;
    v_out = u_out + vB0;
end

Answer 2

Я попытался обрисовать несколько простой подход, так сказать, обратной стороны конверта, разделенный на четыре простых шага.

Предположим, у вас есть начальные позиции и скорости xA0, vA0 и xB0, vB0 космического корабля. А и Б соответственно. Как вы сказали, B движется без ускорения и с постоянной скоростью vB0. Поэтому он равномерно движется по прямой. Его движение описывается следующим образом: xB = xB0 + t*vB0 Космический корабль A может включать и выключать ускорение постоянной величины a0, но может менять направление по своему усмотрению.

Я действительно надеюсь, что ваш предел скорости удовлетворяет norm(vA0 - vB0) < v_max в противном случае управление ускорением, которое вы должны построить, становится более сложным.

Шаг 1: Убейте разницу между скоростями A и B. Примените постоянное ускорение

a = a0 *(vB0 - vA0) / norm(vB0 - vA0)

к космическому кораблю A. Затем, позиции и скорости A и B изменяются со временем следующим образом:

xA = xA0 + t*vA0 + t^2*a0*(vB0 - vA0)/(2*norm(vB0 - vA0))
vA = vA0 + t*a0*(vB0 - vA0)/norm(vB0 - vA0)
xB = xB0 + t*vB0
vB = vB0

В момент времени t1 = norm(vB0 - vA0)/a0 скорость космического корабля A равна vB0, который по величине и направлению равен скорости космического корабля B. В t1 если A отключает ускорение и удерживает его, он будет двигаться параллельно B, просто со смещением в пространстве.

Объяснение: (не требуется для алгоритма, но объясняет вычисления, использованные в следующих шагах) Поскольку космический корабль B движется равномерно по прямой линии с постоянной скоростью vB0, он фактически определяет инерциальная система координат. Другими словами, если мы переведем исходную систему координат и прикрепим ее к B, новая система движется с постоянной скоростью вдоль прямой линии и поэтому также является инерционной. Преобразование является галилеевым, поэтому можно определить следующее изменение координат (в обоих направлениях)

y = x - xB0 - t*vB0
u = v - vB0

x = y + xB0 + t*vB0
v = u + vB0

Во время t1, начиная с шага 1, позиции двух космических кораблей

xA1 = xA0 + t1*vA0 + t1^2*a0*(vB0 - vA0)/(2*norm(vB0 - vA0))
xB1 = xB0 + t*vB0

и их скорости vA1 = vB1 = vB0. Таким образом,

yA1 = xA1 - xB0 - t1*vB0  

yB1 = xB1 - xB0 - t1*vB0 = xB0 + t1*vB0 - xB0 - t1*vB0  = 0

В этой системе координат, если в момент времени t1 A отключает свое ускорение и удерживает его, он будет просто неподвижен, то есть его положение yA1 не изменится со временем. Теперь все, что нам нужно сделать, это переместить А из точки yA1 в 0 вдоль отрезка прямой AB, определяемого вектором - yA1 = vector(AB) (указывающего от А к началу В). Идея состоит в том, что теперь А может просто двигаться с постоянным ускорением по AB в течение некоторого времени (t2-t1), набирая некоторую скорость uA2, которая не превышает ваш предел скорости morm(uA2 + vB0) < v_max, затем отключить ускорение и летать в течение некоторого периода времени (t3-t2), которое должно быть определено, со скоростью uA2 и, наконец, включается замедление по AB для времени (t4-t3) = (t2-t1), а во время t4 A и B встречаются, и скорость A равна 0 (в новой системе координат, летающей с B). Это означает, что два корабля находятся в одном месте и имеют одинаковую скорость (как вектор) в исходной системе координат.

Теперь

yA = yA1 - (t-t1)^2*a0*yA1/(2*norm(yA1))
uA = (t-t1)*a0*yA1/norm(yA1)

, поэтому на t2 (все точки yA1, yA2, yA3 и 0 коллинеарны):

yA2 = yA1 - (t2-t1)^2*a0*yA1/(2*norm(yA1)) = (norm(yA1)-(t2-t1)^2*a0/(2*norm(yA1))) * yA1
uA2 = (t2-t1)*a0*yA1/norm(yA1)

norm(yA2 - yA1) = norm( yA1 - (t2-t1)^2*a0*yA1/(2*norm(yA1)) - yA1 ) 
                = norm(- (t2-t1)^2*a0*yA1/(2*norm(yA1))) 
                = (t2-t1)^2*(a0/2)*norm(yA1/norm(yA1))
                = (t2-t1)^2*a0/2
norm(yA1) = norm(yA2 - yA1) + norm(yA3 - yA2) + norm(0 - yA3)

norm(yA3 - yA2) = norm(yA1) - norm(yA2 - yA1) - norm(0 - yA3) 
                =  norm(yA1) - (t2-t1)^2*a0

(t3-t2) = norm(yA3 - yA2) / norm(uA2) = ( norm(yA1) - (t2-t1)^2*a0 )/norm(uA2)

Теперь вернемся к исходной системе координат.

yA1 = xA1 - xB1
uA2 = vA2 - vB0 
(t3-t2) = ( norm(xA1 - xB1) - (t2-t1)^2*a0 )/norm(vA2 - vB0)

, поэтому важный расчет здесь таков: как только вы выберете t2, вы получите для вычисления

t3 = t2 + ( norm(xA1 - xB1) - (t2-t1)^2*a0 )/norm(vA2 - vB0)

Шаг 2 : Как уже упоминалось, в момент времени t1 с шага 1 позиции двух космических кораблей равны

xA1 = xA0 + t1*vA0 + t1^2*a0*(vB0 - vA0)/(2*norm(vB0 - vA0))
xB1 = xB0 + t*vB0

, а их скорости равны vA1 = vB1 = vB0.

At время t1 применить ускорение a = a0*(xB1 - xA1)/norm(xB1 - xA1). Затем положения и скорости A и B меняются со временем следующим образом:

xA = xA1 + (t-t1)*vB0 + (t-t1)^2*a0*(xB1 - xA1)/(2*norm(xB1 - xA1))
vA = vB0 + (t-t1)*a0*(xB1 - xA1)/norm(xB1 - xA1)
xB = xB1 + (t-t1)*vB0 or if you prefer xB = xB0 + t*vB0
vB = vB0

Выберите любое t2, которое удовлетворяет

t2 <= t1 + sqrt( norm(xA1 - xB1)/a0 )   (the time to get to the middle of ``AB`` accelerating)

and such that it satisfies

norm( vB0 - (t2 - t1)*a0*(xA1 - xB1)/norm(xA1 - xB1) ) < v_max

Тогда в момент t2 вы получите позиции скорости

xA2 = xA1 + (t2-t1)*vB0 + (t2-t1)^2*a0*(xB1 - xA1)/(2*norm(xB1 - xA1))
vA2 = vB0 + (t2-t1)*a0*(xB1 - xA1)/norm(xB1 - xA1)
xB2 = xB1 + (t2-t1)*vB0   or if you prefer xB2 = xB0 + t2*vB0
vB2 = vB0

Шаг 3: Рассчитать следующий момент времени

t3 = t2 + ( norm(xA1 - xB1) - (t2-t1)^2*a0 )/norm(vA2 - vB0)

и, поскольку А движется с постоянной скоростью vA2 по прямой линия:

xA3 = xA2 + (t3-t2)*vA2
vA3 = vA2
xB3 = xB2 + (t3-t2)*vB0   or if you prefer xB3 = xB0 + t3*vB0
vB3 = vB0

Шаг 4: Это последний отрезок, когда А замедляется, чтобы встретиться с В:

t4 = t3 + (t2-t1)

В момент t3 применить ускорение a = a0*(xA1 - xB1)/norm(xA1 - XB1), прямо противоположный шагу 2. Затем положения и скорости A и B меняются со временем следующим образом:

xA = xA3 + (t-t3)*vB3 + (t-t3)^2*a0*(xA1 - xB1)/(2*norm(xA1 - xB1))
vA = vB3 + (t-t3)*a0*(xA1 - xB1)/norm(xA1 - xB1)
xB = xB3 + (t-t3)*vB0 or if you prefer xB = xB0 + t*vB0
vB = vB0

и для t4 мы должны иметь

xA4 = xB4  and vA4 = vB0

Теперь я понимаю, что существует справедливая количество деталей, так что возможно у меня есть некоторые опечатки и, возможно, ошибки. Тем не менее, идея кажется мне правильной, но я советую вам повторить некоторые расчеты, чтобы быть уверенным.

Answer 3

Предположим, у вас есть начальные позиции и скорости xA0, vA0 и xB0, vB0 космических кораблей A и B соответственно. Как вы сказали, B движется без ускорения и с постоянной скоростью vB0. Поэтому он равномерно движется по прямой. Его движение описывается как: xB = xB0 + t*vB0. Космический корабль A может включать и выключать ускорение постоянной величины a0, но может менять направление по своему усмотрению.

Поскольку космический корабль B движется равномерно по прямой линии с постоянной скоростью vB0, он фактически определяет инерциальную систему координат. Другими словами, если мы переведем исходную систему координат и прикрепим ее к B, новая система движется с постоянной скоростью вдоль прямой линии и поэтому также является инерционной. Преобразование является галилеевым, поэтому можно определить следующее изменение координат (в обоих направлениях)

y = x - xB0 - t*vB0
u = v - vB0

x = y + xB0 + t*vB0
v = u + vB0

, в частности, для B для любого момента времени t мы получим

yB = xB - xB0 - t*vB0 = xB0 + t*vB0 - xB0 - t*vB0  = 0``

В момент t=0,

yA0 = xA0 - xB0  

uA0 = vA0 - vB0

Итак, мы собираемся спроектировать элемент управления в этой новой системе координат, и они переместят его обратно в исходную. Во-первых, космический корабль А движется так, что его скорость всегда имеет одинаковую величину norm(uA) = norm(uA0), но его направление меняется равномерно. Для этого можно просто взять вектор перекрестных произведений

L0 = cross(yA0, uA0) / ( norm( cross(yA0, uA0) ) * norm(uA0) )

и в каждый момент времени t применить ускорение

a = a0 * cross(L0, uA)

Это означает, что закон движения А удовлетворяет дифференциальным уравнениям

dyA/dt = uA

duA/dt = a0 * cross(L0 , uA)

, затем

d/dt (dot(uA, uA)) = 2 * dot(uA, duA/dt) = 2 * dot(uA, a0 * cross(L0 , uA))
                  = 2 * a0 * dot(uA, cross(L0 , uA)) 
                  = 0

, что возможно только тогда, когда norm(uA)^2 = dot(uA, uA) = norm(uA0), т. Е. Величина norm(uA) = norm(uA0) для всех t постоянна.

Давайте проверим норму величины ускорения:

norm(a) = a0 * norm( cross(L0, uA)) = a0 * norm(L0) * norm(uA)
        = a0 * norm( cross(yA0, uA0)/( norm( cross(yA0, uA0) )*norm(uA0) ) )*norm(uA0)
        = a0 * norm( cross(yA0, uA0) )/( norm( cross(yA0, uA0) )*norm(uA0) ) )*norm(uA0) 
        = a0

Так как norm(uA) = norm(uA0) = const вершина скорости A, нарисованная как вектор uA от начала координат B, всегда лежит на сфере norm(uA) = norm(uA0) с центром в начале координат. В то же время

d/dt ( dot(L0, uA) ) = dot(L0, duA/dt) = a0 * dot(L0, cross(L0, uA)) = 0
which means that
dot(L0, uA) = const = dot(L0, uA0) = 0

следовательно uA всегда лежит на плоскости, перпендикулярной вектору L0 и проходящей через начало координат. Таким образом, uA указывает на пересечение указанной плоскости со сферой norm(uA) = norm(uA0), то есть uA пересекает круг. Другими словами, уравнение

duA/dt = a0 cross(L0, uA) 

определяет поворот вокруг начала вектора uA в плоскости, проходящей через начало координат и перпендикулярной L0. Затем возьмите

dyA/dt - a0*cross(L0, yA) = uA - a0*cross(L0, yA) 

и дифференцируйте его по t:

duA/dt - a0*cross(L0, dyA/dt) = duA/dt - a0*cross(L0, uA) = 0 

, что означает, что существует постоянный вектор, такой что dyA/dt - a0*cross(L0, yA) = const_vect, и мы можем переписать это последнее уравнение как

dyA/dt = a0*cross(L0, yA - cA)

and even like

d/dt( yA - cA ) = a0*cross(L0, yA - cA)

, что по тем же аргументам, что и для uA, означает, что yA - cA пересекает окружность с центром в начале координат и в плоскости, перпендикулярной L0. Следовательно, yA пересекает круг в плоскости через начало координат, перпендикулярно L0 и центрировано в cA. Нужно просто найти радиус и центр круга. Тогда движение А по уравнениям

dyA/dt = uA

duA/dt = a0* cross(L0, uA)

сводится к уравнению

dyA/dt = a0 * cross(L0, yA - cA)

yA(0) = yA0

Чтобы найти радиус R, мы устанавливаем время t=0:

uA0 = a0 * cross(L0, yA0 - cA)
so
norm(uA0) = a0 * norm(cross(L0, yA0 - cA)) = a0 * norm(L0) * norm(yA0 - cA)
          = a0 * norm(L0) * R
norm(L0) = 1 / norm(uA0)

R = norm(uA0)^2 / a0

Центр расположен вдоль вектора, перпендикулярного как uA0, так и L0, поэтому

cA = yA0 + R * cross(L0, uA0) / (norm(L0)*norm(uA0))

Затем мы можем установить 2D систему координат в плоскости, в которой движение происходит путем выбора начала координат yA0 и единичных перпендикулярных векторов uA0/norm(uA0) и -cross(L0, uA0) / (norm(L0)*norm(uA0)). Таким образом, движение A в системе координат, движущееся равномерно по прямой линии с B, можно описать как

yA(t) = yA0  + R * sin(a0 * t / norm(L0)) * uA0 / norm(uA0)
             - R * cos(a0 * t / norm(L0)) * cross(L0, uA0) / (norm(L0)*norm(uA0))

, что является решением проблемы начальных значений:

dyA/dt = uA0

duA/dt = a0 * cross(L0, uA)

yA(0) = yA0
uA(0) = uA0

Так Мое новое предложение состоит в том, чтобы включить

Шаг 0: В течение периода времени t от 0 до t0 применить следующее ускорение и движение, которое вращает направление скорости вектор A:

yA0 = xA0 - xB0
uA0 = vA0 - vB0
L0 = cross(yA0, uA0) / ( norm( cross(yA0, uA0) ) * norm(uA0) )
a = a0 * cross(L0, uA0)
R = norm(uA0)^2 / a0

yA(t) = yA0 + R * cos(a0*t/norm(uA0)) / (norm(L0)*norm(uA0))
            + R * sin(a0*t/norm(uA0)) * uA0/norm(uA0)
            - R * cos(a0*t/norm(uA0)) * cross(L0, uA0) / (norm(L0)*norm(uA0))

xA(t) = yA(t) + xB0 + t * vB0 = 
      = xA0 + t * vB0 + R * cos(a0*t/norm(uA0)) / (norm(L0)*norm(uA0))
            + R * sin(a0*t/norm(uA0)) * uA0/norm(uA0)
            - R * cos(a0*t/norm(uA0)) * cross(L0, uA0) / (norm(L0)*norm(uA0))

до момента времени t0, выбранный так, чтобы направление скорости vA(t) находилось в лучшем положении относительно vB0, так что с момента t0 вы Можно применить четыре шага, изложенные в моем предыдущем ответе. Конечно, вы также можете использовать это новое управление круговыми движениями, чтобы создать свою собственную комбинацию, которая вам больше нравится.