diff --git a/benchmark/power_cone_benchmarks.py b/benchmark/power_cone_benchmarks.py
new file mode 100644
index 0000000000..cd9d587ef1
--- /dev/null
+++ b/benchmark/power_cone_benchmarks.py
@@ -0,0 +1,652 @@
+"""
+Copyright, the CVXPY authors
+Licensed under the Apache License, Version 2.0 (the "License");
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at
+    https://www.apache.org/licenses/LICENSE-2.0
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License.
+
+Benchmarks for power cone atoms comparing approx=True vs approx=False.
+
+- approx=True: Uses SOC approximation via rational decomposition
+- approx=False: Uses native power cones (requires solver support)
+
+These benchmarks measure compilation and solve times for:
+- power(x, p) with various exponents
+- geo_mean(x, p) with various weight configurations
+- pnorm(x, p) with various p values
+"""
+
+import cvxpy as cp
+import numpy as np
+
+
+# =============================================================================
+# Power Atom Benchmarks
+# =============================================================================
+
+class PowerApproxTrue:
+    """Power atom with approx=True (SOC approximation)."""
+
+    def setup(self):
+        np.random.seed(42)
+        n = 5000
+        self.x = cp.Variable(n, pos=True)
+        c = np.random.rand(n) + 0.1
+
+        # Power with p=1.5 (convex, non-integer)
+        objective = cp.Minimize(cp.sum(c @ self.x))
+        constraints = [
+            cp.sum(cp.power(self.x, 1.5, approx=True)) <= 100,
+            self.x >= 0.01,
+        ]
+        self.problem = cp.Problem(objective, constraints)
+
+    def time_compile_problem(self):
+        self.problem.get_problem_data(solver=cp.CLARABEL)
+
+    def time_solve_problem(self):
+        self.problem.solve(solver=cp.CLARABEL, verbose=False)
+
+
+class PowerApproxFalse:
+    """Power atom with approx=False (native power cones)."""
+
+    def setup(self):
+        np.random.seed(42)
+        n = 5000
+        self.x = cp.Variable(n, pos=True)
+        c = np.random.rand(n) + 0.1
+
+        # Power with p=1.5 (convex, non-integer)
+        objective = cp.Minimize(cp.sum(c @ self.x))
+        constraints = [
+            cp.sum(cp.power(self.x, 1.5, approx=False)) <= 100,
+            self.x >= 0.01,
+        ]
+        self.problem = cp.Problem(objective, constraints)
+
+    def time_compile_problem(self):
+        self.problem.get_problem_data(solver=cp.CLARABEL)
+
+    def time_solve_problem(self):
+        self.problem.solve(solver=cp.CLARABEL, verbose=False)
+
+
+class PowerNegativeExponentApproxTrue:
+    """Power atom with negative exponent, approx=True."""
+
+    def setup(self):
+        np.random.seed(42)
+        n = 1000
+        self.x = cp.Variable(n, pos=True)
+        c = np.random.rand(n) + 0.1
+
+        # Power with p=-0.5 (convex, negative exponent)
+        objective = cp.Minimize(cp.sum(cp.power(self.x, -0.5, approx=True)))
+        constraints = [
+            c @ self.x >= 10,
+            self.x >= 0.01,
+            self.x <= 10,
+        ]
+        self.problem = cp.Problem(objective, constraints)
+
+    def time_compile_problem(self):
+        self.problem.get_problem_data(solver=cp.CLARABEL)
+
+    def time_solve_problem(self):
+        self.problem.solve(solver=cp.CLARABEL, verbose=False)
+
+
+class PowerNegativeExponentApproxFalse:
+    """Power atom with negative exponent, approx=False."""
+
+    def setup(self):
+        np.random.seed(42)
+        n = 1000
+        self.x = cp.Variable(n, pos=True)
+        c = np.random.rand(n) + 0.1
+
+        # Power with p=-0.5 (convex, negative exponent)
+        objective = cp.Minimize(cp.sum(cp.power(self.x, -0.5, approx=False)))
+        constraints = [
+            c @ self.x >= 10,
+            self.x >= 0.01,
+            self.x <= 10,
+        ]
+        self.problem = cp.Problem(objective, constraints)
+
+    def time_compile_problem(self):
+        self.problem.get_problem_data(solver=cp.CLARABEL)
+
+    def time_solve_problem(self):
+        self.problem.solve(solver=cp.CLARABEL, verbose=False)
+
+
+class PowerFractionalExponentApproxTrue:
+    """Power atom with fractional exponent (0 < p < 1), approx=True."""
+
+    def setup(self):
+        np.random.seed(42)
+        n = 1000
+        self.x = cp.Variable(n, pos=True)
+        c = np.random.rand(n) + 0.1
+
+        # Power with p=0.3 (concave)
+        objective = cp.Maximize(cp.sum(cp.power(self.x, 0.3, approx=True)))
+        constraints = [
+            c @ self.x <= 100,
+            self.x >= 0.01,
+        ]
+        self.problem = cp.Problem(objective, constraints)
+
+    def time_compile_problem(self):
+        self.problem.get_problem_data(solver=cp.CLARABEL)
+
+    def time_solve_problem(self):
+        self.problem.solve(solver=cp.CLARABEL, verbose=False)
+
+
+class PowerFractionalExponentApproxFalse:
+    """Power atom with fractional exponent (0 < p < 1), approx=False."""
+
+    def setup(self):
+        np.random.seed(42)
+        n = 1000
+        self.x = cp.Variable(n, pos=True)
+        c = np.random.rand(n) + 0.1
+
+        # Power with p=0.3 (concave)
+        objective = cp.Maximize(cp.sum(cp.power(self.x, 0.3, approx=False)))
+        constraints = [
+            c @ self.x <= 100,
+            self.x >= 0.01,
+        ]
+        self.problem = cp.Problem(objective, constraints)
+
+    def time_compile_problem(self):
+        self.problem.get_problem_data(solver=cp.CLARABEL)
+
+    def time_solve_problem(self):
+        self.problem.solve(solver=cp.CLARABEL, verbose=False)
+
+
+# =============================================================================
+# Geometric Mean Benchmarks
+# =============================================================================
+
+class GeoMeanApproxTrue:
+    """Geometric mean with approx=True (SOC approximation)."""
+
+    def setup(self):
+        np.random.seed(42)
+        n = 100
+        m = 50
+        self.x = cp.Variable(n, pos=True)
+        A = np.random.rand(m, n) + 0.1
+        b = np.random.rand(m) * 10
+
+        objective = cp.Maximize(cp.geo_mean(self.x, approx=True))
+        constraints = [
+            A @ self.x <= b,
+            self.x >= 0.01,
+        ]
+        self.problem = cp.Problem(objective, constraints)
+
+    def time_compile_problem(self):
+        self.problem.get_problem_data(solver=cp.CLARABEL)
+
+    def time_solve_problem(self):
+        self.problem.solve(solver=cp.CLARABEL, verbose=False)
+
+
+class GeoMeanApproxFalse:
+    """Geometric mean with approx=False (native power cones)."""
+
+    def setup(self):
+        np.random.seed(42)
+        n = 100
+        m = 50
+        self.x = cp.Variable(n, pos=True)
+        A = np.random.rand(m, n) + 0.1
+        b = np.random.rand(m) * 10
+
+        objective = cp.Maximize(cp.geo_mean(self.x, approx=False))
+        constraints = [
+            A @ self.x <= b,
+            self.x >= 0.01,
+        ]
+        self.problem = cp.Problem(objective, constraints)
+
+    def time_compile_problem(self):
+        self.problem.get_problem_data(solver=cp.CLARABEL)
+
+    def time_solve_problem(self):
+        self.problem.solve(solver=cp.CLARABEL, verbose=False)
+
+
+class GeoMeanLargeApproxTrue:
+    """Large geometric mean with approx=True."""
+
+    def setup(self):
+        np.random.seed(42)
+        n = 500
+        m = 250
+        self.x = cp.Variable(n, pos=True)
+        A = np.random.rand(m, n) + 0.1
+        b = np.random.rand(m) * 10
+
+        objective = cp.Maximize(cp.geo_mean(self.x, approx=True))
+        constraints = [
+            A @ self.x <= b,
+            self.x >= 0.01,
+        ]
+        self.problem = cp.Problem(objective, constraints)
+
+    def time_compile_problem(self):
+        self.problem.get_problem_data(solver=cp.CLARABEL)
+
+    def time_solve_problem(self):
+        self.problem.solve(solver=cp.CLARABEL, verbose=False)
+
+
+class GeoMeanLargeApproxFalse:
+    """Large geometric mean with approx=False."""
+
+    def setup(self):
+        np.random.seed(42)
+        n = 500
+        m = 250
+        self.x = cp.Variable(n, pos=True)
+        A = np.random.rand(m, n) + 0.1
+        b = np.random.rand(m) * 10
+
+        objective = cp.Maximize(cp.geo_mean(self.x, approx=False))
+        constraints = [
+            A @ self.x <= b,
+            self.x >= 0.01,
+        ]
+        self.problem = cp.Problem(objective, constraints)
+
+    def time_compile_problem(self):
+        self.problem.get_problem_data(solver=cp.CLARABEL)
+
+    def time_solve_problem(self):
+        self.problem.solve(solver=cp.CLARABEL, verbose=False)
+
+
+class GeoMeanWeightedApproxTrue:
+    """Weighted geometric mean with approx=True."""
+
+    def setup(self):
+        np.random.seed(42)
+        n = 100
+        m = 50
+        self.x = cp.Variable(n, pos=True)
+        A = np.random.rand(m, n) + 0.1
+        b = np.random.rand(m) * 10
+
+        # Non-uniform positive weights
+        weights = np.random.rand(n) + 0.5
+        weights = weights / weights.sum()
+
+        objective = cp.Maximize(cp.geo_mean(self.x, p=weights, approx=True))
+        constraints = [
+            A @ self.x <= b,
+            self.x >= 0.01,
+        ]
+        self.problem = cp.Problem(objective, constraints)
+
+    def time_compile_problem(self):
+        self.problem.get_problem_data(solver=cp.CLARABEL)
+
+    def time_solve_problem(self):
+        self.problem.solve(solver=cp.CLARABEL, verbose=False)
+
+
+class GeoMeanWeightedApproxFalse:
+    """Weighted geometric mean with approx=False."""
+
+    def setup(self):
+        np.random.seed(42)
+        n = 100
+        m = 50
+        self.x = cp.Variable(n, pos=True)
+        A = np.random.rand(m, n) + 0.1
+        b = np.random.rand(m) * 10
+
+        # Non-uniform positive weights
+        weights = np.random.rand(n) + 0.5
+        weights = weights / weights.sum()
+
+        objective = cp.Maximize(cp.geo_mean(self.x, p=weights, approx=False))
+        constraints = [
+            A @ self.x <= b,
+            self.x >= 0.01,
+        ]
+        self.problem = cp.Problem(objective, constraints)
+
+    def time_compile_problem(self):
+        self.problem.get_problem_data(solver=cp.CLARABEL)
+
+    def time_solve_problem(self):
+        self.problem.solve(solver=cp.CLARABEL, verbose=False)
+
+
+# =============================================================================
+# P-norm Benchmarks
+# =============================================================================
+
+class PnormApproxTrue:
+    """P-norm with p=1.5, approx=True."""
+
+    def setup(self):
+        np.random.seed(42)
+        n = 5000
+        m = 2500
+        self.x = cp.Variable(n)
+        A = np.random.randn(m, n)
+        b = np.random.randn(m)
+
+        objective = cp.Minimize(cp.pnorm(self.x, p=1.5, approx=True))
+        constraints = [
+            A @ self.x == b,
+        ]
+        self.problem = cp.Problem(objective, constraints)
+
+    def time_compile_problem(self):
+        self.problem.get_problem_data(solver=cp.CLARABEL)
+
+    def time_solve_problem(self):
+        self.problem.solve(solver=cp.CLARABEL, verbose=False)
+
+
+class PnormApproxFalse:
+    """P-norm with p=1.5, approx=False."""
+
+    def setup(self):
+        np.random.seed(42)
+        n = 5000
+        m = 2500
+        self.x = cp.Variable(n)
+        A = np.random.randn(m, n)
+        b = np.random.randn(m)
+
+        objective = cp.Minimize(cp.pnorm(self.x, p=1.5, approx=False))
+        constraints = [
+            A @ self.x == b,
+        ]
+        self.problem = cp.Problem(objective, constraints)
+
+    def time_compile_problem(self):
+        self.problem.get_problem_data(solver=cp.CLARABEL)
+
+    def time_solve_problem(self):
+        self.problem.solve(solver=cp.CLARABEL, verbose=False)
+
+
+class PnormP3ApproxTrue:
+    """P-norm with p=3, approx=True."""
+
+    def setup(self):
+        np.random.seed(42)
+        n = 5000
+        m = 2500
+        self.x = cp.Variable(n)
+        A = np.random.randn(m, n)
+        b = np.random.randn(m)
+
+        objective = cp.Minimize(cp.pnorm(self.x, p=3, approx=True))
+        constraints = [
+            A @ self.x == b,
+        ]
+        self.problem = cp.Problem(objective, constraints)
+
+    def time_compile_problem(self):
+        self.problem.get_problem_data(solver=cp.CLARABEL)
+
+    def time_solve_problem(self):
+        self.problem.solve(solver=cp.CLARABEL, verbose=False)
+
+
+class PnormP3ApproxFalse:
+    """P-norm with p=3, approx=False."""
+
+    def setup(self):
+        np.random.seed(42)
+        n = 5000
+        m = 2500
+        self.x = cp.Variable(n)
+        A = np.random.randn(m, n)
+        b = np.random.randn(m)
+
+        objective = cp.Minimize(cp.pnorm(self.x, p=3, approx=False))
+        constraints = [
+            A @ self.x == b,
+        ]
+        self.problem = cp.Problem(objective, constraints)
+
+    def time_compile_problem(self):
+        self.problem.get_problem_data(solver=cp.CLARABEL)
+
+    def time_solve_problem(self):
+        self.problem.solve(solver=cp.CLARABEL, verbose=False)
+
+
+class PnormFractionalApproxTrue:
+    """P-norm with p=0.5, approx=True."""
+
+    def setup(self):
+        np.random.seed(42)
+        n = 1000
+        m = 500
+        self.x = cp.Variable(n, pos=True)
+        A = np.random.rand(m, n) + 0.1
+        b = np.random.rand(m) * 10
+        c = np.random.rand(n)
+
+        # Minimize linear, constrain pnorm
+        objective = cp.Minimize(c @ self.x)
+        constraints = [
+            cp.pnorm(self.x, p=0.5, approx=True) >= 10,
+            A @ self.x <= b,
+            self.x >= 0.01,
+        ]
+        self.problem = cp.Problem(objective, constraints)
+
+    def time_compile_problem(self):
+        self.problem.get_problem_data(solver=cp.CLARABEL)
+
+    def time_solve_problem(self):
+        self.problem.solve(solver=cp.CLARABEL, verbose=False)
+
+
+class PnormFractionalApproxFalse:
+    """P-norm with p=0.5, approx=False."""
+
+    def setup(self):
+        np.random.seed(42)
+        n = 1000
+        m = 500
+        self.x = cp.Variable(n, pos=True)
+        A = np.random.rand(m, n) + 0.1
+        b = np.random.rand(m) * 10
+        c = np.random.rand(n)
+
+        # Minimize linear, constrain pnorm
+        objective = cp.Minimize(c @ self.x)
+        constraints = [
+            cp.pnorm(self.x, p=0.5, approx=False) >= 10,
+            A @ self.x <= b,
+            self.x >= 0.01,
+        ]
+        self.problem = cp.Problem(objective, constraints)
+
+    def time_compile_problem(self):
+        self.problem.get_problem_data(solver=cp.CLARABEL)
+
+    def time_solve_problem(self):
+        self.problem.solve(solver=cp.CLARABEL, verbose=False)
+
+
+# =============================================================================
+# SCS Solver Benchmarks
+# =============================================================================
+
+class PowerSCSApproxTrue:
+    """Power atom with SCS solver, approx=True."""
+
+    def setup(self):
+        np.random.seed(42)
+        n = 1000
+        self.x = cp.Variable(n, pos=True)
+        c = np.random.rand(n) + 0.1
+
+        objective = cp.Minimize(cp.sum(c @ self.x))
+        constraints = [
+            cp.sum(cp.power(self.x, 1.5, approx=True)) <= 100,
+            self.x >= 0.01,
+        ]
+        self.problem = cp.Problem(objective, constraints)
+
+    def time_compile_problem(self):
+        self.problem.get_problem_data(solver=cp.SCS)
+
+    def time_solve_problem(self):
+        self.problem.solve(solver=cp.SCS, verbose=False)
+
+
+class PowerSCSApproxFalse:
+    """Power atom with SCS solver, approx=False."""
+
+    def setup(self):
+        np.random.seed(42)
+        n = 1000
+        self.x = cp.Variable(n, pos=True)
+        c = np.random.rand(n) + 0.1
+
+        objective = cp.Minimize(cp.sum(c @ self.x))
+        constraints = [
+            cp.sum(cp.power(self.x, 1.5, approx=False)) <= 100,
+            self.x >= 0.01,
+        ]
+        self.problem = cp.Problem(objective, constraints)
+
+    def time_compile_problem(self):
+        self.problem.get_problem_data(solver=cp.SCS)
+
+    def time_solve_problem(self):
+        self.problem.solve(solver=cp.SCS, verbose=False)
+
+
+class GeoMeanSCSApproxTrue:
+    """Geometric mean with SCS solver, approx=True."""
+
+    def setup(self):
+        np.random.seed(42)
+        n = 100
+        m = 50
+        self.x = cp.Variable(n, pos=True)
+        A = np.random.rand(m, n) + 0.1
+        b = np.random.rand(m) * 10
+
+        objective = cp.Maximize(cp.geo_mean(self.x, approx=True))
+        constraints = [
+            A @ self.x <= b,
+            self.x >= 0.01,
+        ]
+        self.problem = cp.Problem(objective, constraints)
+
+    def time_compile_problem(self):
+        self.problem.get_problem_data(solver=cp.SCS)
+
+    def time_solve_problem(self):
+        self.problem.solve(solver=cp.SCS, verbose=False)
+
+
+class GeoMeanSCSApproxFalse:
+    """Geometric mean with SCS solver, approx=False."""
+
+    def setup(self):
+        np.random.seed(42)
+        n = 100
+        m = 50
+        self.x = cp.Variable(n, pos=True)
+        A = np.random.rand(m, n) + 0.1
+        b = np.random.rand(m) * 10
+
+        objective = cp.Maximize(cp.geo_mean(self.x, approx=False))
+        constraints = [
+            A @ self.x <= b,
+            self.x >= 0.01,
+        ]
+        self.problem = cp.Problem(objective, constraints)
+
+    def time_compile_problem(self):
+        self.problem.get_problem_data(solver=cp.SCS)
+
+    def time_solve_problem(self):
+        self.problem.solve(solver=cp.SCS, verbose=False)
+
+
+if __name__ == '__main__':
+    import time
+
+    benchmarks = [
+        # Power atom benchmarks
+        ("PowerApproxTrue", PowerApproxTrue),
+        ("PowerApproxFalse", PowerApproxFalse),
+        ("PowerNegativeExponentApproxTrue", PowerNegativeExponentApproxTrue),
+        ("PowerNegativeExponentApproxFalse", PowerNegativeExponentApproxFalse),
+        ("PowerFractionalExponentApproxTrue", PowerFractionalExponentApproxTrue),
+        ("PowerFractionalExponentApproxFalse", PowerFractionalExponentApproxFalse),
+        # Geometric mean benchmarks
+        ("GeoMeanApproxTrue", GeoMeanApproxTrue),
+        ("GeoMeanApproxFalse", GeoMeanApproxFalse),
+        ("GeoMeanLargeApproxTrue", GeoMeanLargeApproxTrue),
+        ("GeoMeanLargeApproxFalse", GeoMeanLargeApproxFalse),
+        ("GeoMeanWeightedApproxTrue", GeoMeanWeightedApproxTrue),
+        ("GeoMeanWeightedApproxFalse", GeoMeanWeightedApproxFalse),
+        # P-norm benchmarks
+        ("PnormApproxTrue", PnormApproxTrue),
+        ("PnormApproxFalse", PnormApproxFalse),
+        ("PnormP3ApproxTrue", PnormP3ApproxTrue),
+        ("PnormP3ApproxFalse", PnormP3ApproxFalse),
+        ("PnormFractionalApproxTrue", PnormFractionalApproxTrue),
+        ("PnormFractionalApproxFalse", PnormFractionalApproxFalse),
+        # SCS solver benchmarks
+        ("PowerSCSApproxTrue", PowerSCSApproxTrue),
+        ("PowerSCSApproxFalse", PowerSCSApproxFalse),
+        ("GeoMeanSCSApproxTrue", GeoMeanSCSApproxTrue),
+        ("GeoMeanSCSApproxFalse", GeoMeanSCSApproxFalse),
+    ]
+
+    print("Running power cone benchmarks...")
+    print("Comparing approx=True (SOC) vs approx=False (power cones)\n")
+    print(f"{'Benchmark':<40} {'Compile (s)':<15} {'Solve (s)':<15}")
+    print("=" * 70)
+
+    for name, cls in benchmarks:
+        try:
+            bench = cls()
+            bench.setup()
+
+            # Time compilation
+            start = time.time()
+            bench.time_compile_problem()
+            compile_time = time.time() - start
+
+            # Time solve
+            start = time.time()
+            bench.time_solve_problem()
+            solve_time = time.time() - start
+
+            print(f"{name:<40} {compile_time:<15.4f} {solve_time:<15.4f}")
+        except Exception as e:
+            print(f"{name:<40} ERROR: {e}")
+
+    print("\nDone!")
diff --git a/benchmark/power_cone_iterations.py b/benchmark/power_cone_iterations.py
new file mode 100644
index 0000000000..b612b00158
--- /dev/null
+++ b/benchmark/power_cone_iterations.py
@@ -0,0 +1,208 @@
+"""
+Power cone iteration count experiments.
+
+Compares solve performance between native power cones (approx=False)
+and SOC approximation (approx=True) for various atoms.
+
+Usage:
+    uv run python power_cone_iterations.py CLARABEL
+    uv run python power_cone_iterations.py SCS
+    uv run python power_cone_iterations.py MOSEK
+"""
+
+import argparse
+import cvxpy as cp
+import numpy as np
+
+
+def get_problem_stats(prob, solver):
+    """Extract problem statistics after solving."""
+    data = prob.get_problem_data(solver)
+    cones = data[0]['dims']
+    cone_info = []
+
+    # ConeDims attributes: zero, nonneg, exp, soc, psd, p3d, pnd
+    # SCS also has 'p' as an alias for p3d
+    if hasattr(cones, 'zero') and cones.zero > 0:
+        cone_info.append(f"Zero: {cones.zero}")
+    if hasattr(cones, 'nonneg') and cones.nonneg > 0:
+        cone_info.append(f"Nonneg: {cones.nonneg}")
+    if hasattr(cones, 'exp') and cones.exp > 0:
+        cone_info.append(f"Exp: {cones.exp}")
+    if hasattr(cones, 'soc') and cones.soc:
+        cone_info.append(f"SOC: {len(cones.soc)} cones")
+    if hasattr(cones, 'psd') and cones.psd:
+        cone_info.append(f"PSD: {len(cones.psd)} cones")
+    if hasattr(cones, 'p3d') and cones.p3d:
+        cone_info.append(f"Pow3D: {len(cones.p3d)} cones")
+    if hasattr(cones, 'pnd') and cones.pnd:
+        cone_info.append(f"GenPow: {len(cones.pnd)} cones")
+    # SCS uses 'p' for power cones (list of alpha values)
+    if hasattr(cones, 'p') and cones.p:
+        cone_info.append(f"Pow3D: {len(cones.p)} cones")
+
+    n_vars = data[0]['A'].shape[1]
+    n_constrs = data[0]['A'].shape[0]
+
+    return n_vars, n_constrs, ", ".join(cone_info)
+
+
+def print_results(label, prob, solver, solver_name):
+    """Print formatted results after solving."""
+    n_vars, n_constrs, cone_info = get_problem_stats(prob, solver)
+
+    # Get solver stats - all solvers have num_iters and solve_time
+    iterations = prob.solver_stats.num_iters
+    solve_time = prob.solver_stats.solve_time
+
+    # Handle None iterations (some solvers don't report this)
+    if iterations is None:
+        iterations = 0
+        time_per_iter = 0
+    else:
+        time_per_iter = (solve_time / iterations * 1000) if iterations > 0 else 0
+
+    print(f"\n{label}:")
+    print(f"  Variables:      {n_vars}")
+    print(f"  Constraints:    {n_constrs}")
+    print(f"  Cones:          {cone_info}")
+    print(f"  Status:         {prob.status}")
+    print(f"  Iterations:     {iterations}")
+    print(f"  Solve time:     {solve_time*1000:.2f} ms")
+    print(f"  Time/iteration: {time_per_iter:.2f} ms")
+    print(f"  Optimal value:  {prob.value:.6f}")
+
+
+def run_geo_mean_experiment(solver, solver_name, n=200):
+    """Compare geo_mean with power cones vs SOC approximation."""
+    print(f"\n{'='*60}")
+    print(f"geo_mean n={n} with {solver_name}")
+    print('='*60)
+
+    np.random.seed(42)
+
+    # Power cones (approx=False)
+    x = cp.Variable(n, pos=True)
+    obj = cp.sum(x)
+    constraint = [cp.geo_mean(x, approx=False) >= 0.1, x <= 10]
+    prob = cp.Problem(cp.Minimize(obj), constraint)
+    prob.solve(solver=solver)
+    print_results("approx=False (Power Cones)", prob, solver, solver_name)
+
+    # SOC approximation (approx=True)
+    x2 = cp.Variable(n, pos=True)
+    obj2 = cp.sum(x2)
+    constraint2 = [cp.geo_mean(x2, approx=True) >= 0.1, x2 <= 10]
+    prob2 = cp.Problem(cp.Minimize(obj2), constraint2)
+    prob2.solve(solver=solver)
+    print_results("approx=True (SOC)", prob2, solver, solver_name)
+
+
+def run_pnorm_experiment(solver, solver_name, n=2000, p=1.5):
+    """Compare pnorm with power cones vs SOC approximation."""
+    print(f"\n{'='*60}")
+    print(f"pnorm(x, {p}) n={n} with {solver_name}")
+    print('='*60)
+
+    np.random.seed(42)
+    c = np.random.randn(n)
+
+    # Power cones (approx=False)
+    x = cp.Variable(n)
+    obj = c @ x
+    constraint = [cp.pnorm(x, p, approx=False) <= 1]
+    prob = cp.Problem(cp.Minimize(obj), constraint)
+    prob.solve(solver=solver)
+    print_results("approx=False (Power Cones)", prob, solver, solver_name)
+
+    # SOC approximation (approx=True)
+    x2 = cp.Variable(n)
+    obj2 = c @ x2
+    constraint2 = [cp.pnorm(x2, p, approx=True) <= 1]
+    prob2 = cp.Problem(cp.Minimize(obj2), constraint2)
+    prob2.solve(solver=solver)
+    print_results("approx=True (SOC)", prob2, solver, solver_name)
+
+
+def run_power_experiment(solver, solver_name, n=1000, p=0.5):
+    """Compare power with power cones vs SOC approximation."""
+    print(f"\n{'='*60}")
+    print(f"power(x, {p}) n={n} with {solver_name}")
+    print('='*60)
+
+    np.random.seed(42)
+    c = np.random.randn(n)
+
+    # Power cones (approx=False)
+    x = cp.Variable(n, pos=True)
+    obj = c @ x
+    constraint = [cp.sum(cp.power(x, p, approx=False)) >= n, x <= 10]
+    prob = cp.Problem(cp.Minimize(obj), constraint)
+    prob.solve(solver=solver)
+    print_results("approx=False (Power Cones)", prob, solver, solver_name)
+
+    # SOC approximation (approx=True)
+    x2 = cp.Variable(n, pos=True)
+    obj2 = c @ x2
+    constraint2 = [cp.sum(cp.power(x2, p, approx=True)) >= n, x2 <= 10]
+    prob2 = cp.Problem(cp.Minimize(obj2), constraint2)
+    prob2.solve(solver=solver)
+    print_results("approx=True (SOC)", prob2, solver, solver_name)
+
+
+def run_inv_prod_experiment(solver, solver_name, n=50):
+    """Compare inv_prod with power cones vs SOC approximation."""
+    print(f"\n{'='*60}")
+    print(f"inv_prod n={n} with {solver_name}")
+    print('='*60)
+
+    np.random.seed(42)
+    c = np.random.randn(n)
+
+    # Power cones (approx=False)
+    x = cp.Variable(n, pos=True)
+    obj = c @ x
+    constraint = [cp.inv_prod(x, approx=False) <= 1, x >= 0.1, x <= 10]
+    prob = cp.Problem(cp.Minimize(obj), constraint)
+    prob.solve(solver=solver)
+    print_results("approx=False (Power Cones)", prob, solver, solver_name)
+
+    # SOC approximation (approx=True)
+    x2 = cp.Variable(n, pos=True)
+    obj2 = c @ x2
+    constraint2 = [cp.inv_prod(x2, approx=True) <= 1, x2 >= 0.1, x2 <= 10]
+    prob2 = cp.Problem(cp.Minimize(obj2), constraint2)
+    prob2.solve(solver=solver)
+    print_results("approx=True (SOC)", prob2, solver, solver_name)
+
+
+def main():
+    parser = argparse.ArgumentParser(description="Power cone iteration experiments")
+    parser.add_argument("solver", choices=["CLARABEL", "SCS", "MOSEK"],
+                        help="Solver to use")
+    parser.add_argument("--experiment", "-e",
+                        choices=["geo_mean", "pnorm", "power", "inv_prod", "all"],
+                        default="all", help="Which experiment to run")
+    parser.add_argument("--n", type=int, default=None,
+                        help="Problem size (default varies by experiment)")
+    args = parser.parse_args()
+
+    solver = getattr(cp, args.solver)
+    solver_name = args.solver
+
+    experiments = {
+        "geo_mean": lambda: run_geo_mean_experiment(solver, solver_name, n=args.n or 200),
+        "pnorm": lambda: run_pnorm_experiment(solver, solver_name, n=args.n or 2000),
+        "power": lambda: run_power_experiment(solver, solver_name, n=args.n or 1000),
+        "inv_prod": lambda: run_inv_prod_experiment(solver, solver_name, n=args.n or 50),
+    }
+
+    if args.experiment == "all":
+        for exp in experiments.values():
+            exp()
+    else:
+        experiments[args.experiment]()
+
+
+if __name__ == "__main__":
+    main()