다중 객체 추적기용 추적 필터 자동으로 조정하기
이 예제에서는 trackingFilterTuner 객체를 사용하여 추적 필터를 자동으로 조정하는 방법을 보여줍니다. 조정 후에는 다중 타깃 추적기에서 조정 결과를 사용하여 추적기의 추적 성능을 개선합니다.
조정되지 않은 필터로 추적기 실행하기
이 예제는 Air Traffic Control 예제의 시나리오를 기반으로 합니다.
필터 조정이 필요한지 여부를 조사하기 위해 디폴트 필터 초기화 함수를 사용하여 trackerGNN System object를 생성합니다. 또한 OSPA(Optimal Sub-Pattern Assignment: 최적의 부분패턴 연관) 메트릭 trackOSPAMetric을 정의하여 추적 결과를 정량적으로 분석합니다. 참고로, OSPA 값이 낮을수록 추적 품질이 더 좋습니다. 할당 임계값을 조정하려면 슬라이더를 사용합니다.
tracker = trackerGNN(AssignmentThreshold =100); metric = trackOSPAMetric(Distance = "posabserr", CutoffDistance = 200);
원래 예제와 비교했을 때, 레이다는 빔이 더 넓고 레이다 스캔 속도가 더 빠르도록 수정되어 시나리오의 각 타깃에 대해 더 많은 검출을 생성할 수 있습니다.
load("ATCScenario.mat","scenario"); helperRunScenario(scenario, tracker, metric);
![Figure contains an axes object. The axes object with title OSPA Metric, xlabel Time [seconds], ylabel OSPA contains an object of type line.](../../examples/fusion/win64/AutomaticallyTuneATrackingFilterExample_02.png){kind=tracking}
결과를 보면, 할당 임계값이 가장 높게 설정된 경우에도 추적기의 성능은 기대에 미치지 못하며 모든 타깃을 추적할 수는 없습니다.
시나리오로부터 조정 데이터 만들기
필터를 자동으로 조정하려면 truth 객체의 이력과 각 truth 객체별 검출 로그와 함께 trackingFilterTuner 객체를 사용합니다. truth 데이터는 timetable 객체로 구성된 셀형 배열로, 각 셀은 하나의 truth 객체의 내역을 나타냅니다. 검출 로그는 셀형 배열입니다. 검출 로그의 번째 셀은 truth 데이터의 번째 truth 객체에 대응하는 objectDetection 객체로 구성된 셀형 배열입니다.
시나리오로부터 조정 데이터를 만들려면 monteCarloRun 함수를 사용합니다. 이 함수는 여러 다른 잡음 구현으로 시나리오 기록을 만듭니다. 기록을 조정 데이터로 변환하려면 예제에 첨부된 helperRecordingToTuningData 지원 함수를 사용합니다. 간결한 설명을 위해 이 예제에서는 시나리오에 대해 몬테카를로 구현을 2회 실행하여 만들어진 조정 데이터가 포함된 파일을 직접 불러옵니다.
% recording = monteCarloRun(scenario,2); % [detections, truthTable] = recordingToTuningData(recording); load("ATCTuningData.mat","detections","truthTable"); helperVisualizeTuningData(detections,truthTable);
디폴트 조정 파라미터로 필터 조정하기
추적기가 빠르게 움직이는 플랫폼을 추적하기 어려운 이유를 파악하려면 추적기에서 사용되는 디폴트 initcvekf 함수로 초기화되는 필터를 살펴봅니다.
sampleDetection = detections{1}{1};
filter = initcvekf(sampleDetection);필터는 등속도 모션 모델 규칙을 사용하며 상태 벡터는 로 정의됩니다. 여기서 , , 는 사각형 프레임의 위치 성분이고 , , 는 속도 성분입니다. 필터의 공정 잡음은 객체의 속도 변화(즉, 가속도)에 대한 불확실성을 나타냅니다.
disp(filter.ProcessNoise);
1 0 0
0 1 0
0 0 1
마찬가지로, 상태 공분산은 위치 성분과 속도 성분의 불확실성을 나타냅니다. initcvekf 함수는 검출 측정 잡음을 직접 사용하여 상태 불확실성 공분산의 위치 성분을 초기화합니다. 그러나 속도 성분은 측정에 보고되지 않으므로 이러한 속도 성분 각각의 불확실성은 기본적으로 100으로 설정되며 이는 10m/s의 표준편차를 나타냅니다. 여객기의 경우 불확실성이 당연히 더 높을 것입니다.
disp(filter.StateCovariance);
1.0e+04 *
0.9582 0 -0.1154 0 -0.0824 0
0 0.0100 0 0 0 0
-0.1154 0 0.0488 0 -0.3422 0
0 0 0 0.0100 0 0
-0.0824 0 -0.3422 0 4.1085 0
0 0 0 0 0 0.0100
trackingFilterTuner 객체를 만듭니다.
tuner = trackingFilterTuner; disp(tuner)
trackingFilterTuner with properties:
FilterInitializationFcn: "initcvekf"
TunablePropertiesSource: "Default"
Cost: "RMSE"
UseMex: 0
UseParallel: 0
Solver: "active-set"
Display: "Text"
기본적으로, initcvekf 함수로 초기화되는 trackingEKF 객체를 조정할 때 조정기는 필터의 ProcessNoise 속성을 조정하고 필터 추정값과 truth 간의 RMSE(RMS 오차)를 비용으로 사용합니다. 일반적으로 공정 잡음이 정의하기가 가장 어렵습니다. 그 이유는 객체 모션이 필터에 사용되는 모션 모델과 어떻게 다른지에 대한 정보가 거의 없는 경우가 많기 때문입니다. 그다음으로 조정하기 어려운 것은 초기 상태 공분산이며, 이는 다음 섹션에서 조정할 것입니다.
조정을 가속화하려면 체크박스를 사용하여 UseMex 속성을 true로 설정합니다. 이 설정을 위해서는 MATLAB Coder 라이선스가 필요합니다. 이 라이선스를 사용할 수 없는 경우 체크박스를 선택 취소하십시오.
tuner.UseMex = 
true;마찬가지로, 병렬 처리를 사용하려면 UseParallel 속성을 true로 설정합니다. 이 경우 Parallel Computing Toolbox 라이선스가 필요합니다. 이 라이선스를 사용할 수 없는 경우 체크박스를 선택 취소하십시오.
tuner.UseParallel = 
true;fmincon 최적화 알고리즘, particleswarm 최적화 알고리즘, patternsearch 최적화 알고리즘을 기반으로 하는 3개의 솔버 중 하나를 사용할 수 있습니다. 최적화 문제의 복잡도, 실행에 허용된 시간, 기타 기준에 따라 최적화 알고리즘을 선택하십시오. Optimization Toolbox™와 Global Optimization Toolbox™ 문서를 참조하여 어떤 솔버를 사용할지 결정할 수 있습니다.
'fmincon' 솔버를 사용하려면 Optimization Toolbox 라이선스가 필요합니다.
'particleswarm' 솔버 또는 'patternsearch' 솔버를 사용하려면 Global Optimization Toolbox 라이선스가 필요합니다.
solver = 
"fmincon";시간 경과에 따라 해결 진행 상황을 반복적으로 표시하도록 솔버 옵션을 설정합니다. 슬라이더를 사용하여 솔버가 사용할 수 있는 최대 반복 횟수를 제어합니다. 반복 횟수를 높게 설정하면 조정기에 더 많은 시간이 필요합니다. 모든 타깃을 추적하도록 필터를 조정하려면 모든 타깃에 대해 조정 데이터를 사용합니다.
tuner.Solver = solver; tuner.SolverOptions = optimoptions(solver, Display = "iter", ... MaxIter =15); tic; [tunedParams, bestCost] = tune(tuner,detections,truthTable);
Starting parallel pool (parpool) using the 'Processes' profile ...
Connected to parallel pool with 6 workers.
Generating code.
Code generation successful.
Your initial point x0 is not between bounds lb and ub; FMINCON
shifted x0 to strictly satisfy the bounds.
Iter RMSE Step Size
0 78.4801 0.0000
First-order Norm of
Iter F-count f(x) Feasibility optimality step
0 7 7.848006e+01 0.000e+00 1.613e-01
1 14 7.839574e+01 0.000e+00 1.852e-01 3.538e-01
1 78.3957 0.3538
2 21 7.810211e+01 0.000e+00 1.639e-01 1.463e+00
2 78.1021 1.4628
3 28 7.768567e+01 0.000e+00 1.813e-01 3.044e+00
3 77.6857 3.0440
4 35 7.720462e+01 0.000e+00 2.038e-01 4.543e+00
4 77.2046 4.5434
5 43 7.703294e+01 0.000e+00 2.229e-01 1.866e+00
5 77.0329 1.8661
6 50 7.686042e+01 0.000e+00 2.013e-01 1.277e+00
6 76.8604 1.2772
7 57 7.654589e+01 0.000e+00 1.457e-01 1.737e+00
7 76.5459 1.7370
8 64 7.603169e+01 0.000e+00 6.552e-02 4.124e+00
8 76.0317 4.1242
9 71 7.600035e+01 0.000e+00 6.463e-02 8.611e-01
9 76.0003 0.8611
10 78 7.589348e+01 0.000e+00 6.237e-02 1.662e+00
10 75.8935 1.6618
11 85 7.592252e+01 0.000e+00 4.813e-02 5.195e-01
11 75.9225 0.5195
12 92 7.591004e+01 0.000e+00 4.924e-02 1.346e-01
12 75.9100 0.1346
13 99 7.588873e+01 0.000e+00 3.311e-02 3.816e-01
13 75.8887 0.3816
14 106 7.586347e+01 0.000e+00 3.171e-02 1.262e+00
14 75.8635 1.2622
15 113 7.584819e+01 0.000e+00 2.707e-02 1.731e+00
15 75.8482 1.7310
Solver stopped prematurely.
fmincon stopped because it exceeded the iteration limit,
options.MaxIterations = 1.500000e+01.
disp("Time it took to tune the filter: " + toc);Time it took to tune the filter: 135.5945
disp("Optimized cost: " + bestCost);Optimized cost: 75.8482
최적화 과정이 비용을 크게 개선하지 못했다는 것을 알 수 있습니다. 이 경우에는 높은 비용의 주요 동인이 부적절한 초기 속도 공분산(조정에 포함되지 않았던 항목)이기 때문에 예상된 일입니다.
trackingEKF 객체에 대해 조정된 속성은 2개의 필드 ProcessNoise와 StateCovariance를 갖는 구조체로 반환됩니다. 해당 값은 아래에 표시되어 있습니다. 상태 공분산은 조정되지 않았습니다.
disp(tunedParams.ProcessNoise);
111.6332 28.3525 96.7912 28.3525 48.7698 13.3867 96.7912 13.3867 97.0984
disp(tunedParams.StateCovariance);
1.0e+04 *
0.9582 0 -0.1154 0 -0.0824 0
0 0.0100 0 0 0 0
-0.1154 0 0.0488 0 -0.3422 0
0 0 0 0.0100 0 0
-0.0824 0 -0.3422 0 4.1085 0
0 0 0 0 0 0.0100
조정된 필터를 가져오려면 필터의 각 속성을 구조체에서 그에 대응되는 조정된 값으로 설정하거나 필터의 setTunedProperties 객체 함수를 사용합니다.
filter = feval(tuner.FilterInitializationFcn, sampleDetection); setTunedProperties(filter,tunedParams); disp(filter.ProcessNoise);
111.6332 28.3525 96.7912 28.3525 48.7698 13.3867 96.7912 13.3867 97.0984
plotFilterErrors 객체 함수를 사용하여 추정 정확도를 표시합니다. 플롯에는 모든 모션 차원에서의 위치 추정 오차와 속도 추정 오차가 표시됩니다. 각 선은 필터를 한 번 실행했을 때의 추정값과, 연결된 truth 간의 오차를 나타냅니다. 대역은 필터 상태 공분산으로부터 얻은 해당 차원에서의 3표준편차("3시그마")를 나타냅니다. 편향되지 않은 필터의 경우, 필터 추정값 오차의 평균이 0에 가까워야 합니다.
또한 일관된 결과를 내는 가우스 필터의 추정값은 약 99%의 경우 3시그마 범위 내에 있어야 합니다. 필터가 과다신뢰하는 경우 3표준편차 범위의 위반은 더 많아질 것입니다. 필터가 과소신뢰하는 경우 실제 필터 오차는 3표준편차 범위에 비해 매우 작을 것입니다. 추적에서는 불확실성에 대한 올바른 척도를 제공하는 일관된 결과를 내는 필터를 갖추는 것이 중요합니다.
plotFilterErrors(tuner);
![Figure Filter estimate error contains 6 axes objects. Axes object 1 with xlabel Time [sec], ylabel X_{Error} contains 12 objects of type line, patch. Axes object 2 with xlabel Time [sec], ylabel V_X_{Error} contains 12 objects of type line, patch. Axes object 3 with xlabel Time [sec], ylabel Y_{Error} contains 12 objects of type line, patch. Axes object 4 with xlabel Time [sec], ylabel V_Y_{Error} contains 12 objects of type line, patch. Axes object 5 with xlabel Time [sec], ylabel Z_{Error} contains 12 objects of type line, patch. Axes object 6 with xlabel Time [sec], ylabel V_Z_{Error} contains 12 objects of type line, patch. These objects represent Realizations, 3\sigma.](../../examples/fusion/win64/AutomaticallyTuneATrackingFilterExample_04.png){kind=tracking}
조정을 사용자 지정하기
공정 잡음만 조정하는 경우 필터 추정값 오차는 3표준편차 범위를 훨씬 초과합니다. 특히 X축과 Y축 방향의 속도에 대해서는 더욱 그렇습니다. 그 이유는 속도에 대한 초기 상태 공분산 값이 너무 낮기 때문입니다. 다시 말해 필터가 속도 추정값을 너무 과대평가한다는 뜻입니다. 또한 조정 알고리즘이 큰 초기 오차를 보정하려고 시도하기 때문에 조정 후에 공정 잡음 값이 너무 높아집니다. 그 결과 추정값은 필요 이상으로 측정값을 더 가깝게 따라가고 필터 오차는 필터가 수렴된 후에도 큽니다.
이 문제를 해결하기 위해 사용자 지정 속성을 사용하도록 조정기를 설정합니다. 필터 조정 가능형 속성을 만들려면 추적에 사용되는 것과 동일한 모션 모델로 같은 유형의 추적 필터 객체를 만듭니다. 그런 다음 필터의 tunableProperties 객체 함수를 사용하여 tunableFilterProperties 객체를 생성합니다.
filter = feval(tuner.FilterInitializationFcn, sampleDetection); tfp = tunableProperties(filter); disp(tfp)
Tunable properties for object of type: trackingEKF
Property: ProcessNoise
PropertyValue: [1 0 0;0 1 0;0 0 1]
TunedQuantity: Square root
IsTuned: true
TunedQuantityValue: [1 0 0;0 1 0;0 0 1]
TunableElements: [1 4 5 7 8 9]
LowerBound: [0 0 0 0 0 0]
UpperBound: [10 10 10 10 10 10]
Property: StateCovariance
PropertyValue: [9582.02019480753 0 -1154.41787165046 0 -823.540813608376 0;0 100 0 0 0 0;-1154.41787165046 0 488.321270547992 0 -3422.28418831314 0;0 0 0 100 0 0;-823.540813608376 0 -3422.28418831314 0 41085.3267410573 0;0 0 0 0 0 100]
TunedQuantity: Square root of initial value
IsTuned: false
필터 조정 가능형 속성을 사용하여, 조정해야 할 속성과 각 속성을 조정하는 방법을 정의합니다. IsTuned 플래그를 사용하여 속성의 조정 가능성을 설정합니다. 이 경우 StateCovariance 속성이 조정되도록 설정합니다.
마찬가지로, 조정할 요소를 정의합니다. 6×6 상태 공분산 행렬이 양의 준정부호 대칭 행렬이 되도록 하고 최적화할 파라미터의 개수를 36개에서 21개로 줄이기 위해, 촐레스키 분해를 사용하여 상부 삼각 행렬과 그 행렬의 전치의 곱으로 공분산을 인코딩합니다. 그래도 여전히 요소의 수가 매우 많기 때문에 솔버에서 상당한 노력이 필요할 것입니다. 한편, 필터 초기화 함수는 이미 검출 측정 잡음을 사용하여 초기 상태 공분산 위치 불확실성을 설정하므로 속도 요소만 조정하면 됩니다. 게다가 위의 결과는 Z 성분의 속도 오차가 이미 적절하게 포함되어 있음을 보여줍니다.
등속도 모델에 사용되는 상태 벡터는 로 정의된다는 것을 기억하십시오. 여기서 , , 는 사각형 프레임의 위치 성분이고 , , 는 속도 성분입니다. 따라서 X 방향과 Y 방향을 따르는 속도 성분의 불확실성을 조정하려면 촐레스키 형식의 (2,2) 요소, (2,4) 요소, (4,4) 요소만 조정해야 합니다.
조정기를 추가로 지원하려면 각 요소의 하한과 상한을 정의합니다. 위의 그래프에는 초당 300미터 미만의 초기 속도 오차가 표시되어 있으므로 여기서는 하한과 상한을 각각 0과 300으로 정의합니다. 하한 벡터와 상한 벡터는 조정 가능형 요소 개수와 동일한 개수의 요소를 가져야 합니다. 여기서는 3입니다.
elements = sub2ind([6 6], [2 2 4], [2 4 4]); % Tune just the XY velocity covariance elements lb = zeros(1,numel(elements)); ub = 300*ones(1,numel(elements)); setPropertyTunability(tfp, "StateCovariance", IsTuned = true, ... TunableElements = elements, LowerBound = lb, UpperBound = ub);
마찬가지로, 공정 잡음의 모든 요소를 조정할 이유가 없습니다. 공정 잡음은 3×3 행렬이며, X 공정 성분과 Y 공정 성분을 제어하는 요소는 (1,1), (1,2), (2,2)입니다. 조정할 요소의 개수를 줄이면 조정기가 더 빠르게 실행될 수 있습니다.
여객기는 대부분 직선 수평 비행을 하기 때문에 상한은 10을 사용합니다.
elements = sub2ind([3 3], [1 1 2], [1 2 2]); % Tune just the XY process elements lb = zeros(1,numel(elements)); ub = 10*ones(1,numel(elements)); setPropertyTunability(tfp, "ProcessNoise", "IsTuned", true, ... "TunableElements", elements, "LowerBound", lb, "UpperBound", ub);
마지막으로, 조정기의 사용자 지정 조정 가능형 속성을 조정 가능형 필터 속성으로 설정합니다.
tuner.TunablePropertiesSource = "Custom";
tuner.CustomTunableProperties = tfp;solver ="fmincon"; tuner.Solver = solver; tuner.SolverOptions = optimoptions(solver, Display = "iter", ... MaxIter =
15); tic; [tunedParams, bestCost] = tune(tuner,detections,truthTable);
Generating code.
Code generation successful.
Your initial point x0 is not between bounds lb and ub; FMINCON
shifted x0 to strictly satisfy the bounds.
Iter RMSE Step Size
0 81.1683 0.0000
First-order Norm of
Iter F-count f(x) Feasibility optimality step
0 7 8.116826e+01 0.000e+00 2.749e+00
1 14 7.303610e+01 0.000e+00 1.375e+00 2.809e+00
1 73.0361 2.8094
2 21 6.933060e+01 0.000e+00 7.783e-01 3.039e+00
2 69.3306 3.0395
3 28 6.597938e+01 0.000e+00 3.954e-01 4.808e+00
3 65.9794 4.8078
4 35 6.375682e+01 0.000e+00 2.266e-01 5.858e+00
4 63.7568 5.8576
5 42 6.215164e+01 0.000e+00 1.328e-01 7.779e+00
5 62.1516 7.7786
6 49 6.112787e+01 0.000e+00 1.123e-01 9.377e+00
6 61.1279 9.3772
7 56 6.052563e+01 0.000e+00 1.073e-01 1.052e+01
7 60.5256 10.5211
8 63 6.021176e+01 0.000e+00 1.184e-01 9.795e+00
8 60.2118 9.7949
9 70 6.012935e+01 0.000e+00 1.139e-01 2.429e+00
9 60.1293 2.4290
10 77 5.990467e+01 0.000e+00 6.043e-02 7.234e+00
10 59.9047 7.2341
11 84 5.986403e+01 0.000e+00 1.322e-01 6.622e+00
11 59.8640 6.6223
12 91 5.984704e+01 0.000e+00 6.894e-02 4.802e-01
12 59.8470 0.4802
13 98 5.982958e+01 0.000e+00 6.043e-02 2.277e+00
13 59.8296 2.2768
14 105 5.980530e+01 0.000e+00 2.285e-02 5.648e+00
14 59.8053 5.6476
15 112 5.979078e+01 0.000e+00 2.116e-02 7.667e+00
15 59.7908 7.6673
Solver stopped prematurely.
fmincon stopped because it exceeded the iteration limit,
options.MaxIterations = 1.500000e+01.
disp("Tuning time: " + toc);Tuning time: 74.622
disp("Best cost: " + bestCost);Best cost: 59.7908
setTunedProperties(filter,tunedParams);
disp("Tuned process noise is:");Tuned process noise is:
disp(filter.ProcessNoise);
34.9996 50.4336 0
50.4336 83.5728 0
0 0 1.0000
disp("Tuned initial state covariance is:");Tuned initial state covariance is:
disp(filter.StateCovariance);
1.0e+04 *
0.9582 0 -0.1154 0 -0.0824 0
0 0.6434 0 0.0231 0 0
-0.1154 0 0.0488 0 -0.3422 0
0 0.0231 0 0.3186 0 0
-0.0824 0 -0.3422 0 4.1085 0
0 0 0 0 0 0.0100
plotFilterErrors(tuner);
![Figure Filter estimate error contains 6 axes objects. Axes object 1 with xlabel Time [sec], ylabel X_{Error} contains 12 objects of type line, patch. Axes object 2 with xlabel Time [sec], ylabel V_X_{Error} contains 12 objects of type line, patch. Axes object 3 with xlabel Time [sec], ylabel Y_{Error} contains 12 objects of type line, patch. Axes object 4 with xlabel Time [sec], ylabel V_Y_{Error} contains 12 objects of type line, patch. Axes object 5 with xlabel Time [sec], ylabel Z_{Error} contains 12 objects of type line, patch. Axes object 6 with xlabel Time [sec], ylabel V_Z_{Error} contains 12 objects of type line, patch. These objects represent Realizations, 3\sigma.](../../examples/fusion/win64/AutomaticallyTuneATrackingFilterExample_05.png){kind=tracking}
몇 가지의 변화를 살펴보면 이러한 조정의 결과가 공정 잡음만 조정한 결과보다 훨씬 낫다는 것을 알 수 있습니다. 첫째, 최적의 RMSE 비용의 값이 이전보다 훨씬 낮습니다. 둘째, 시나리오 시간 전체에 걸쳐 필터 오차가 모두 3시그마 대역 내에 있습니다. 마지막으로, 이러한 오차와 대역이 더 좁아집니다. 이는 필터 추정값이 truth 값에 잘 수렴됨을 나타냅니다.
사용자 지정 비용으로 조정하기
조정기에는 2가지 비용 메트릭인 RMSE와 NEES(Normalized Estimation Error Squared: 정규화된 추정오차제곱)가 내장되어 있습니다. RMSE 비용은 평균절대오차를 최소화하기 위해 노력하지만 해당 추정값을 과다신뢰하거나 과소신뢰하는 필터를 생성할 수도 있습니다. NEES 비용은 [1]에서 제안된 비용을 기반으로 일관된 결과를 내는 필터를 제공하기 위해 노력합니다. RMSE와 NEES는 모두 위치의 오차와 속도의 오차(진리표에 제공된 경우)를 고려합니다.
사용자 지정 비용을 사용하면 더 나은 조정이 가능한 경우가 있습니다. 그러한 경우는 다음과 같습니다.
위치와 속도 이외의 다른 데이터 세트가 포함된 진리표를 사용하는 경우. 또한 조정기는 사용자 지정 비용 사용 시 진리표를 검증하지 않습니다.
내장된 모션 모델(
constvelmscjac,constaccjac,constturnjac또는singer)이 아닌 모션 모델을 사용하는 경우.비용에 추가 요소 또는 RMSE 메트릭과 NEES 메트릭에 정의된 메트릭 이외의 다른 메트릭을 포함해야 하는 경우.
이 예제에 첨부된 helperCostNormalized 함수는 Blackman과 Popoli에 의해 제안된 증강 마할라노비스 거리(Mahalanobis distance)를 사용합니다[2]. 비용은 과소신뢰된 필터 추정값(큰 상태 공분산 값)과 과다신뢰된 필터 추정값(작은 상태 공분산 값) 모두에 벌점을 적용하여 일관된 결과를 내는 필터를 생성하기 위해 노력한다는 점에서 NEES와 비슷합니다.
비용 메트릭을 사용자 지정 메트릭으로 변경합니다. 코드를 다시 생성할 필요는 없습니다.
tuner.Cost = "Custom"; tuner.CustomCostFcn = @(trkHistory, truth) helperCostNormalized(trkHistory, truth, "constvel"); tic; [tunedParams, bestCost] = tune(tuner,detections,truthTable);
Your initial point x0 is not between bounds lb and ub; FMINCON
shifted x0 to strictly satisfy the bounds.
Iter Custom Cost Step Size
0 59.7241 0.0000
First-order Norm of
Iter F-count f(x) Feasibility optimality step
0 7 5.972415e+01 0.000e+00 3.938e+00
1 14 4.250375e+01 0.000e+00 1.074e+00 5.781e+00
1 42.5038 5.7813
2 21 3.997240e+01 0.000e+00 7.540e-01 2.021e+00
2 39.9724 2.0206
3 28 3.647123e+01 0.000e+00 3.655e-01 4.797e+00
3 36.4712 4.7973
4 35 3.470021e+01 0.000e+00 2.080e-01 4.774e+00
4 34.7002 4.7739
5 42 3.348664e+01 0.000e+00 1.137e-01 5.990e+00
5 33.4866 5.9903
6 49 3.273653e+01 0.000e+00 9.179e-02 6.968e+00
6 32.7365 6.9678
7 56 3.232928e+01 0.000e+00 1.024e-01 6.683e+00
7 32.3293 6.6835
8 63 3.218352e+01 0.000e+00 1.054e-01 3.375e+00
8 32.1835 3.3747
9 70 3.212378e+01 0.000e+00 1.031e-01 1.088e+00
9 32.1238 1.0879
10 77 3.186992e+01 0.000e+00 9.164e-02 6.811e+00
10 31.8699 6.8106
11 84 3.163440e+01 0.000e+00 2.814e-02 1.208e+01
11 31.6344 12.0845
12 91 3.164332e+01 0.000e+00 2.814e-02 7.510e-01
12 31.6433 0.7510
13 98 3.163372e+01 0.000e+00 2.814e-02 8.293e-01
13 31.6337 0.8293
14 105 3.157808e+01 0.000e+00 4.563e-02 4.516e+00
14 31.5781 4.5164
15 112 3.149118e+01 0.000e+00 1.438e-01 8.802e+00
15 31.4912 8.8016
Solver stopped prematurely.
fmincon stopped because it exceeded the iteration limit,
options.MaxIterations = 1.500000e+01.
disp("The time needed to tune the filter after code generation: " + toc);The time needed to tune the filter after code generation: 58.2451
setTunedProperties(filter,tunedParams);
disp("Tuned process noise is:");Tuned process noise is:
disp(filter.ProcessNoise);
69.5499 79.3440 0
79.3440 96.3635 0
0 0 1.0000
disp("Tuned initial state covariance is:");Tuned initial state covariance is:
disp(filter.StateCovariance);
1.0e+04 *
0.9582 0 -0.1154 0 -0.0824 0
0 0.3330 0 0.0280 0 0
-0.1154 0 0.0488 0 -0.3422 0
0 0.0280 0 0.3875 0 0
-0.0824 0 -0.3422 0 4.1085 0
0 0 0 0 0 0.0100
plotFilterErrors(tuner);
![Figure Filter estimate error contains 6 axes objects. Axes object 1 with xlabel Time [sec], ylabel X_{Error} contains 12 objects of type line, patch. Axes object 2 with xlabel Time [sec], ylabel V_X_{Error} contains 12 objects of type line, patch. Axes object 3 with xlabel Time [sec], ylabel Y_{Error} contains 12 objects of type line, patch. Axes object 4 with xlabel Time [sec], ylabel V_Y_{Error} contains 12 objects of type line, patch. Axes object 5 with xlabel Time [sec], ylabel Z_{Error} contains 12 objects of type line, patch. Axes object 6 with xlabel Time [sec], ylabel V_Z_{Error} contains 12 objects of type line, patch. These objects represent Realizations, 3\sigma.](../../examples/fusion/win64/AutomaticallyTuneATrackingFilterExample_06.png){kind=tracking}
disp("Optimized cost: " + bestCost);Optimized cost: 31.4912
조정된 필터를 추적기에 사용하기
조정된 필터가 제공하는 추적 품질의 개선 사항을 관찰하기 위해 조정 결과를 다중 타깃 추적기에 사용합니다. 먼저, 조정기의 exportToFunction 객체 함수를 사용하여 필터 초기화 함수를 내보냅니다.
exportToFunction(tuner, 'myTunedInitFcn');조정된 초기화 함수로 생성되는 필터 객체의 조정된 속성을 확인합니다.
filter = myTunedInitFcn(sampleDetection); disp(filter.ProcessNoise);
69.5499 79.3440 0
79.3440 96.3635 0
0 0 1.0000
disp(filter.StateCovariance);
1.0e+04 *
0.9582 0 -0.1154 0 -0.0824 0
0 0.3330 0 0.0280 0 0
-0.1154 0 0.0488 0 -0.3422 0
0 0.0280 0 0.3875 0 0
-0.0824 0 -0.3422 0 4.1085 0
0 0 0 0 0 0.0100
이 필터 초기화 함수를 사용하여 추적기를 생성하고 추적 결과를 관찰합니다. 초기 상태 불확실성을 보다 잘 처리하도록 필터를 조정하면 AssignmentThreshold 속성의 값을 줄일 수 있습니다. 할당 임계값이 작을수록 추적기가 트랙에 검출을 더 잘 할당할 수 있습니다.
tracker = trackerGNN(FilterInitializationFcn = "myTunedInitFcn", AssignmentThreshold = 50);
helperRunScenario(scenario, tracker, metric);
![Figure contains an axes object. The axes object with title OSPA Metric, xlabel Time [seconds], ylabel OSPA contains an object of type line.](../../examples/fusion/win64/AutomaticallyTuneATrackingFilterExample_08.png){kind=tracking}
위의 OSPA 메트릭 플롯은 훨씬 개선된 추적 결과를 보여줍니다. 메트릭에서, 트랙이 설정되면 시뮬레이션을 시작한 지 약 5초 후에 평균 절대 위치 오차는 약 20~80미터이며, 이는 레이다 측정 정확도를 감안할 때 매우 우수한 정확도입니다. 이 결과는 이전 섹션에 나와 있는 조정의 결과와도 일치합니다.
요약
이 예제에서는 추적 필터 조정기를 사용하여 추적 필터를 자동으로 조정하는 방법을 알아보았습니다. 그리고 조정된 파라미터를 제어하는 방법과 조정에 대한 비용을 선택하는 방법을 알아보았으며 다양한 최적화 알고리즘을 사용해 보았습니다. 또한 조정된 필터를 필터 초기화 함수에 내보내는 방법을 알아보았습니다. 이 함수에서 반환하는 조정된 필터는 다중 타깃 추적기에 사용할 수 있습니다.
참고 문헌
[1] Chen, Z., N. Ahmed, S. Julier, and C. Heckman, “Kalman Filter Tuning with Bayesian Optimization.” ArXiv:1912.08601 [Cs, Eess, Math], Dec. 2019. arXiv.org, http://arxiv.org/abs/1912.08601.
[2] Blackman, S., and R. Popoli. Design and Analysis of Modern Tracking Systems. Artech House Radar Library, Boston, 1999.
지원 함수
이 섹션에 포함되어 있는 함수는 이 예제를 위한 헬퍼 용도이며 향후 릴리스에서 제거되거나 수정되거나 이름이 변경될 수 있습니다.
helperRecordingToTuningData
trackingScenarioRecording 객체로 구성된 배열에서 검출 로그와 truth 데이터를 만듭니다.
function [detlog,truth] = helperRecordingToTuningData(recording) detlog = extractDetections(recording); truth = extractTruth(recording); isemptyDetlog = cellfun(@(dl) isempty(dl), detlog); detlog = detlog(~isemptyDetlog); truth = truth(~isemptyDetlog); end function truth = extractTruth(recording) data = recording.RecordedData; numTruths = numel(data(1).Poses); Time = seconds([data.SimulationTime])'; t = cell(1,numTruths); poses = [data.Poses]; for j = 1:numTruths positions = reshape([poses(j,:).Position],3,[])'; velocities = reshape([poses(j,:).Velocity],3,[])'; t{j} = timetable(Time,positions,velocities,'VariableNames',{'Position','Velocity'}); end truth = repmat(t,1,numel(recording)); end function detlog = extractDetections(recording) numRuns = numel(recording); numTruths = numel(recording(1).RecordedData(1).Poses); detlog = repmat({},1,numTruths * numRuns); logIdx = 0; for runIdx = 1:numRuns data = recording(runIdx).RecordedData; detections = vertcat(data.Detections); for truthIdx = 1:numTruths platID = data(1).Poses(truthIdx).PlatformID; fromThisPlat = cellfun(@(d) d.ObjectAttributes{1}.TargetIndex == platID, detections); logIdx = logIdx + 1; detlog{logIdx} = detections(fromThisPlat); end end end
helperVisualizeTuningData
이 함수는 조정 데이터를 시각화하고 theaterPlot 시각화 객체를 반환합니다.
function helperVisualizeTuningData(detections, truthTable) thp = theaterPlot(AxesUnits=["km","km","km"],XLimits=[-35000 35000],YLimits=[-50000 25000]); trp = trajectoryPlotter(thp,DisplayName="Truth Trajectories"); numTruths = numel(truthTable); pos = cell(1,numTruths); for i = 1:numTruths pos{i} = truthTable{i}.Position; end plotTrajectory(trp,pos); knownPlatforms = []; existingDetPlotters = {}; colors = colororder; for j = 1:numTruths d = vertcat(detections{j}{:}); if ~isempty(d) platIndex = d(1).ObjectAttributes{1}.TargetIndex; if ~any(platIndex == knownPlatforms) dtp = detectionPlotter(thp,DisplayName=("Platform"+platIndex+"Detections"),MarkerEdgeColor=colors(platIndex,:)); knownPlatforms(end+1) = platIndex; %#ok<*AGROW> existingDetPlotters{end+1} = dtp; else dtp = existingDetPlotters{platIndex == knownPlatforms}; end detpos = arrayfun(@(d) d.Measurement, d, UniformOutput=false); plotDetection(dtp,[detpos{:}]'); end end end
helperRunScenario
레코딩을 실행하는 함수로, 객체를 추적하고 결과를 시각화하고, 추적 품질을 분석합니다.
function helperRunScenario(scenario, tracker, metric) % Initialize plotters persistent thp plp dep cvp tap if isempty(thp) thp = theaterPlot(AxesUnits=["km","km","km"],XLimits=[-35000 35000],YLimits=[-50000 25000]); end if isempty(plp) plp = platformPlotter(thp,DisplayName="Platforms"); end if isempty(dep) dep = detectionPlotter(thp,DisplayName="Detections"); end if isempty(cvp) cvp = coveragePlotter(thp,DisplayName="Radar Coverage"); end if isempty(tap) tap = trackPlotter(thp, DisplayName="Tracks", ConnectHistory="on", ColorizeHistory="on"); end % Reset all objects clearPlotterData(thp); restart(scenario); reset(tracker); reset(metric); % Modify the radar to rotate faster and provide more detections radar = scenario.Platforms{1}.Sensors{1}; release(radar); radar.FieldOfView(1) = 14; radar.MaxMechanicalScanRate(1) = 750; % Initialize variables. There are 3215 timestamps in this scenario ospa = zeros(1,3215); trackerTimes = zeros(1,3215); detBuffer = {}; tracks = objectTrack.empty; index = 0; % For repeatable results, use a constant RNG seed r = rng(0,'twister'); h = onCleanup(@() rng(r)); while advance(scenario) time = scenario.SimulationTime; poses = platformPoses(scenario); covcon = coverageConfig(scenario); [detections, sencon] = detect(scenario); if ~sencon.IsScanDone detBuffer = vertcat(detBuffer,detections); else tracks = tracker(detBuffer, time); index = index + 1; ospa(index) = metric(tracks,poses(2:4)); trackerTimes(index) = time; detBuffer = {}; end if isempty(tracks) trkpos = zeros(0,3); trkIDs = string.empty; else trkpos = getTrackPositions(tracks,"constvel"); trkIDs = string([tracks.TrackID]); end if isempty(detBuffer) detpos = zeros(0,3); else detpos = cellfun(@(db) db.Measurement, detBuffer, 'UniformOutput', false); detpos = [detpos{:}]'; end plotPlatform(plp, reshape([poses.Position],3,[])'); plotCoverage(cvp, covcon); plotDetection(dep, detpos); plotTrack(tap, trkpos, trkIDs); end figure; plot(trackerTimes(1:index),ospa(1:index)); ylim([0 200]); title('OSPA Metric'); xlabel('Time [seconds]'); ylabel('OSPA'); end
helperCostNormalized
Blackman과 Popoli에 의해 제안된 증강 마할라노비스 거리를 계산합니다.
시간 스텝 k에서의 상태와 몬테카를로 실행 i 간의 오차 및 해당 시간에서의 상태 추정값이 주어지면(즉, 이고 상태 공분산은 P(k,i)인 경우) Blackman 및 Popoli 비용은 다음과 같습니다.
T는 시간 스텝의 개수입니다. N은 몬테카를로 실행 횟수입니다.
항은 큰 값의 상태 공분산에 벌점을 적용하기 위해 추가되는 반면, 마할라노비스 비용 항은 과다신뢰된 필터 추정값에 벌점을 적용하는 역할을 합니다.
전체 비용은 모든 시간 스텝과 실행에서 의 평균입니다.
function cost = helperCostNormalized(trkHistory, truthTable, model) % Process truthTable to position and velocity truths posTruth = [truthTable.Position]; hasVelocity = ismember('Velocity', truthTable.Properties.VariableNames); if hasVelocity velTruth = [truthTable.Velocity]; end [numSteps, numRuns] = size(trkHistory); c = zeros(numSteps, numRuns, 'like', posTruth); for run = 1:numRuns [posEst, posCovs] = getTrackPositions(trkHistory(:,run),model); [velEst, velCovs] = getTrackVelocities(trkHistory(:,run),model); if hasVelocity poserror = posTruth - posEst; velerror = velTruth - velEst; stateerror = [poserror,velerror]; else stateerror = posTruth - posEst; end for step = 1:numSteps if hasVelocity P = blkdiag(posCovs(:,:,step), velCovs(:,:,step)); else P = posCovs(:,:,step); end c(step,run) = stateerror(step,:) /P * stateerror(step,:)' + log(det(P)); end end cost = mean(c,"all"); end
