3Dof_model_modified

Original code에서 Main_simulation, GCU, Simparameter 파일을 수정하였고
Compute_cvx_Euler , Compute_cvx_Euler_velocity_zero, Verify_infeasible함수가 추가되었음.

1. 개요

GCU 에서 매 Step 마다 비행 시간을 계산하여 최적의 추력을 얻는 알고리즘을 추가하였음.

Algorithm

1. NED 좌표계에서의 초기 위치와 초기 속도를 얻는다.
   
2. 최적의 비행 시간을 구하기 위해서 NED 기준 D축의 위치와 속도를 나눈 절대값을 final time으로 둔다.
   
3. final time으로 Compute_cvx_Euler 함수와 Bisection method를 이용하여 최적의 비행시간을 구한다.
   
4. 최적의 시간을 최종 시간으로 저장한후 Compute_cvx_Euler를 이용하여 얻은 추력을 Export 한다.
   
5. 최종 위치에 근접 했을 때 속력을 줄이기 위해서 Compute_cvx_Euler_zero를 이용하여 얻은 추력을 Export한다.

Main_Simulation에서 남은 시간이 0 이거나 비행체가 지면 아래에 위치 할 때 while문을 중단한다.

---

Problem. 1

위 식을 Convexification을 하여 계산을 한 것이 아래의 Compute_cvx_Euler function이다.

2. Pseudo code

function : Compute_cvx_Euler
      input : Position,Velocity,Nstep
    output : next step thrust value

initialization : Objective = []
                         Constraint = []

begin function
       Objective   +=   Fuel mass
       Constraint   +=   translation dynamics
       Constarint   +=   ( lower bound on thrust magnitude < thrust < upper bound on thrust magnitude )
       Constarint   +=    initial position,initial velocity, final position, final velocity
       if CVX(objective,constraint) is Solved
            Thrust   =   CVX(objective,constraint)
       else
            Thrust   =   0
       Return Thrust.
end function

위 코드를 실행하기 위해서는 http://cvxr.com/cvx/ 에서 CVX module 설치가 필요함.
Compute_cvx_Euler_zero 는 본 함수에서 제한조건에서 최종위치를 제거.

---

Bisection method(at GCU)

    if Altitude < 2m
       Thrust = Compute_cvx_Euler_Velocity_zero(position,velocity,Nstep)
       Final time = Final time - dt
    else
        Upper bound (nstep) = (z-axis position / z-axis velocity) / dt
        Lower bound (nstep) = 0
        Epsilon = 1   (stopping criterion)
        if Compute cvx_Euler(position,velocity,Upper bound) is infeasible

            while True do
                Lower bound = Upper bound
                Upper bound = Upper bound +Upper bound
                if Compute cvx_Euler(position,velocity,Upper bound) is feasible
                    break
                end if
            end while
        end if

    while Upper bound - Lower bound > Epsilon do
        Temp_N = 0.5 * (Upper bound + Lower bound)
        if Compute_cvx_Euler(position,velocity,Temp_N) is infeasible
           Lower bound = Temp_N
        elseif Compute_cvx_Euler(position,velocity,Temp_N) is feasible
           Upper bound = Temp_N
        end if
    end while
    end if

---

3. 결과

Simulation1의 Parameter 는 참고문헌 [1]를 참고하였음.

Simulation 1

Table. 1 Parameters

Parameter	Value	Units	-
Initial position	[0  500  -500]	[m]	NED
Initial velocity	[50  0  50]	[m/s]	NED
Initial mass	15000	[kg]
Vehicle dry mass	10000	[kg]
Fuel mass	5000	[kg]
Final position	[0  0  0]	[m]	NED
Final velocity	[0  0  0]	[m/s]	NED
Maximum Thrust	250,000	[N]
Minimum Thrust	100,000	[N]
dt	0.1	[s]

Optimal final time: 14.3 sec
Computation time: 52 min

Simulation result

Fig. 1.1 Position

Fig. 1.2 Velocity

Fig. 1.3 Thrust

Fig. 1.4 Thrust norm

3D-animation

Simulation 2

Table. 2 Parameters

Parameter	Value	Units	-
Initial position	[500  100  -700]	[m]	NED
Initial velocity	[50  -50  50]	[m/s]	NED
Initial mass	15000	[kg]
Vehicle dry mass	10000	[kg]
Fuel mass	5000	[kg]
Final position	[0  0  0]	[m]	NED
Final velocity	[0  0  0]	[m/s]	NED
Maximum Thrust	250,000	[N]
Minimum Thrust	100,000	[N]
dt	0.1	[s]

Simulation 1에서 초기위치,초기속도를 변경

Optimal final time: 17.2 sec
Computation time: 74 min

Simulation result

Fig. 2.1 Position

fig. 2.2 Velocity

Fig. 2.3 Thrust

Fig. 2.4 Thrust norm

3D-animation

---

주어진 Convex problem 을 계산하기 위해, NED 좌표계에서 Table1와 Table2를 활용하여 컴퓨터 시뮬레이션을 수행하였다. Fig 1.4와 Fig 2.4는 제한조건에서의 추력제한을 만족하는 것을 보여주며, Fig 1.1과 Fig 2.1는 최종위치까지 도달하는 것을 확인할 수 있다.

4. 참고문헌

Michael Szmuk, Behcet Aclkmese, Andrew W. Berning Jr., “Successive Convexification for Fuel-Optimal Powered Landing With Aerodynamic Drag and Non-Convex Constraints,”American Institute of Aeronautics and Astronautics