diff --git a/qiskit_optimization/algorithms/qrao/magic_rounding.py b/qiskit_optimization/algorithms/qrao/magic_rounding.py
index bc867a80..4fc67ab2 100644
--- a/qiskit_optimization/algorithms/qrao/magic_rounding.py
+++ b/qiskit_optimization/algorithms/qrao/magic_rounding.py
@@ -232,7 +232,7 @@ def _unpack_measurement_outcome(
         """
         output_bits = []
         # iterate in order over decision variables
-        for q, op in var2op.values():
+        for _, (q, op) in sorted(var2op.items()):
             # get the decoding outcome index for the variable
             # corresponding to this Pauli op.
             op_index = self._OP_INDICES[vars_per_qubit][str(op.paulis[0])]
diff --git a/releasenotes/notes/fix-qrao-d0402da6a5c40121.yaml b/releasenotes/notes/fix-qrao-d0402da6a5c40121.yaml
new file mode 100644
index 00000000..ec10897d
--- /dev/null
+++ b/releasenotes/notes/fix-qrao-d0402da6a5c40121.yaml
@@ -0,0 +1,4 @@
+---
+fixes:
+  - |
+    Fixed a bug of :meth:`.MagicRounding.round` that may return a wrong decision variable order.
\ No newline at end of file
diff --git a/test/algorithms/qrao/test_magic_rounding.py b/test/algorithms/qrao/test_magic_rounding.py
index 15caae12..98bdf583 100644
--- a/test/algorithms/qrao/test_magic_rounding.py
+++ b/test/algorithms/qrao/test_magic_rounding.py
@@ -14,6 +14,7 @@
 import unittest
 from test.optimization_test_case import QiskitOptimizationTestCase
 
+import networkx as nx
 import numpy as np
 from qiskit.circuit import QuantumCircuit
 from qiskit.primitives import Sampler
@@ -27,6 +28,7 @@
     RoundingContext,
     RoundingResult,
 )
+from qiskit_optimization.applications import Maxcut
 from qiskit_optimization.problems import QuadraticProgram
 
 
@@ -42,6 +44,15 @@ def setUp(self):
         self.problem.binary_var("z")
         self.problem.minimize(linear={"x": 1, "y": 2, "z": 3})
 
+        # trivial maxcut problem (solution: {1}, {0, 2} with cut val 2)
+        graph = nx.Graph()
+        graph.add_nodes_from([0, 1, 2])
+        graph.add_edges_from([(0, 1), (1, 2)])
+        maxcut = Maxcut(graph)
+        self.maxcut_problem = maxcut.to_quadratic_program()
+        self.maxcut_optimal_value = 2.0
+        self.maxcut_optimal_solution = [0, 1, 0]
+
     def test_magic_rounding_constructor(self):
         """Test constructor"""
         sampler = Sampler(options={"shots": 10000, "seed": 42})
@@ -287,6 +298,19 @@ def test_magic_rounding_round_weighted_3_1_qrac(self):
             [0.2672612419124245, 0.5345224838248487, 0.8017837257372733],
         )
 
+    def test_mapping_magic_rounding_result(self):
+        """Test the mapping of magic rounding result bits to variables"""
+        encoding = QuantumRandomAccessEncoding(max_vars_per_qubit=1)
+        encoding.encode(self.maxcut_problem)
+        circuit = encoding.state_preparation_circuit(self.maxcut_optimal_solution)
+        rounding_context = RoundingContext(encoding=encoding, expectation_values=0, circuit=circuit)
+        sampler = Sampler(options={"shots": 1})
+        magic_rounding = MagicRounding(sampler=sampler)
+        rounding_result = magic_rounding.round(rounding_context)
+        solution = rounding_result.samples[0]
+        self.assertEqual(solution.fval, self.maxcut_optimal_value)
+        np.testing.assert_allclose(solution.x, self.maxcut_optimal_solution)
+
     def test_magic_rounding_exceptions(self):
         """Test exceptions in the MagicRounding class"""
         encoding = QuantumRandomAccessEncoding(max_vars_per_qubit=3)