# Comprehensive Examples
This document provides detailed, real-world examples demonstrating how to use `scirs2-optimize` for various optimization problems.
## Table of Contents
1. [Basic Optimization Examples](#basic-optimization-examples)
2. [Machine Learning Applications](#machine-learning-applications)
3. [Engineering Applications](#engineering-applications)
4. [Financial Applications](#financial-applications)
5. [Scientific Computing](#scientific-computing)
6. [Advanced Use Cases](#advanced-use-cases)
## Basic Optimization Examples
### 1. Simple Quadratic Function
The most basic example - minimizing a simple quadratic function.
```rust
use scirs2_optimize::prelude::*;
use ndarray::Array1;
fn main() -> Result<(), Box<dyn std::error::Error>> {
// Minimize f(x, y) = (x-1)² + (y-2)²
// Analytical minimum: (1, 2) with value 0
let objective = |x: &ArrayView1<f64>| -> f64 {
(x[0] - 1.0).powi(2) + (x[1] - 2.0).powi(2)
};
// Analytical gradient: [2(x-1), 2(y-2)]
let gradient = |x: &ArrayView1<f64>| -> Array1<f64> {
Array1::from_vec(vec![2.0 * (x[0] - 1.0), 2.0 * (x[1] - 2.0)])
};
let x0 = Array1::from_vec(vec![0.0, 0.0]);
let mut options = Options::default();
options.jac = Some(Box::new(gradient));
let result = minimize(objective, &x0, Method::BFGS, Some(options))?;
println!("Solution: [{:.6}, {:.6}]", result.x[0], result.x[1]);
println!("Function value: {:.2e}", result.fun);
println!("Iterations: {}", result.nit);
println!("Success: {}", result.success);
// Verify solution
assert!((result.x[0] - 1.0).abs() < 1e-6);
assert!((result.x[1] - 2.0).abs() < 1e-6);
assert!(result.fun < 1e-10);
Ok(())
}
```
### 2. Rosenbrock Function
The classic Rosenbrock function - a challenging test case for optimization algorithms.
```rust
use scirs2_optimize::prelude::*;
use ndarray::Array1;
// Rosenbrock function: f(x,y) = (a-x)² + b(y-x²)²
// Global minimum at (a,a²) with value 0
// Standard parameters: a=1, b=100
fn rosenbrock(x: &ArrayView1<f64>) -> f64 {
let (a, b) = (1.0, 100.0);
(a - x[0]).powi(2) + b * (x[1] - x[0].powi(2)).powi(2)
}
// Analytical gradient
fn rosenbrock_gradient(x: &ArrayView1<f64>) -> Array1<f64> {
let (a, b) = (1.0, 100.0);
let df_dx = -2.0 * (a - x[0]) - 4.0 * b * x[0] * (x[1] - x[0].powi(2));
let df_dy = 2.0 * b * (x[1] - x[0].powi(2));
Array1::from_vec(vec![df_dx, df_dy])
}
fn main() -> Result<(), Box<dyn std::error::Error>> {
let x0 = Array1::from_vec(vec![-1.2, 1.0]); // Standard starting point
// Compare different algorithms
let algorithms = vec![
("BFGS", Method::BFGS),
("L-BFGS", Method::LBFGS),
("CG", Method::CG),
("Powell", Method::Powell),
];
for (name, method) in algorithms {
let mut options = Options::default();
options.max_iter = 1000;
// Provide gradient for gradient-based methods
if matches!(method, Method::BFGS | Method::LBFGS | Method::CG) {
options.jac = Some(Box::new(rosenbrock_gradient));
}
let start = std::time::Instant::now();
let result = minimize(rosenbrock, &x0, method, Some(options))?;
let duration = start.elapsed();
println!("\n{} Results:", name);
println!(" Solution: [{:.6}, {:.6}]", result.x[0], result.x[1]);
println!(" Function value: {:.2e}", result.fun);
println!(" Iterations: {}", result.nit);
println!(" Function evaluations: {}", result.nfev);
println!(" Time: {:.2?}", duration);
println!(" Success: {}", result.success);
}
Ok(())
}
```
### 3. Constrained Optimization Example
Minimize a function subject to constraints.
```rust
use scirs2_optimize::{constrained::*, prelude::*};
use ndarray::Array1;
fn main() -> Result<(), Box<dyn std::error::Error>> {
// Minimize f(x,y) = (x-2)² + (y-1)²
// Subject to: x + y ≤ 3 (inequality)
// x² + y² = 5 (equality)
let objective = |x: &ArrayView1<f64>| -> f64 {
(x[0] - 2.0).powi(2) + (x[1] - 1.0).powi(2)
};
// Inequality constraint: 3 - x - y ≥ 0
let ineq_constraint = |x: &ArrayView1<f64>| -> f64 {
3.0 - x[0] - x[1]
};
// Equality constraint: x² + y² - 5 = 0
let eq_constraint = |x: &ArrayView1<f64>| -> f64 {
x[0].powi(2) + x[1].powi(2) - 5.0
};
let constraints = vec![
Constraint {
fun: Box::new(ineq_constraint),
constraint_type: ConstraintType::Inequality,
jac: None, // Auto-computed
},
Constraint {
fun: Box::new(eq_constraint),
constraint_type: ConstraintType::Equality,
jac: None,
},
];
// Starting point on the circle x² + y² = 5
let x0 = Array1::from_vec(vec![2.0, 1.0]);
let result = minimize_constrained(
objective,
&x0,
&constraints,
Method::SLSQP,
None
)?;
println!("Constrained minimum:");
println!(" Solution: [{:.6}, {:.6}]", result.x[0], result.x[1]);
println!(" Objective value: {:.6}", result.fun);
// Verify constraints
let ineq_val = ineq_constraint(&result.x.view());
let eq_val = eq_constraint(&result.x.view());
println!(" Inequality constraint (≥0): {:.6}", ineq_val);
println!(" Equality constraint (=0): {:.6}", eq_val);
assert!(ineq_val >= -1e-6, "Inequality constraint violated");
assert!(eq_val.abs() < 1e-6, "Equality constraint violated");
Ok(())
}
```
## Machine Learning Applications
### 1. Logistic Regression from Scratch
Implement logistic regression using stochastic optimization.
```rust
use scirs2_optimize::stochastic::*;
use ndarray::{Array1, Array2, ArrayView1, ArrayView2};
use rand::prelude::*;
struct LogisticRegressionData {
features: Array2<f64>,
labels: Array1<f64>,
n_samples: usize,
n_features: usize,
}
impl LogisticRegressionData {
fn new(features: Array2<f64>, labels: Array1<f64>) -> Self {
let (n_samples, n_features) = features.dim();
Self { features, labels, n_samples, n_features }
}
fn sigmoid(z: f64) -> f64 {
1.0 / (1.0 + (-z.max(-709.0).min(709.0)).exp()) // Clamp to prevent overflow
}
fn predict(&self, params: &ArrayView1<f64>, features: &ArrayView1<f64>) -> f64 {
let z = features.dot(params);
Self::sigmoid(z)
}
fn log_likelihood(&self, params: &ArrayView1<f64>, batch_indices: &[f64]) -> f64 {
let mut ll = 0.0;
for &idx in batch_indices {
let i = (idx as usize) % self.n_samples;
let features = self.features.row(i);
let y = self.labels[i];
let p = self.predict(params, &features);
let p_clamped = p.max(1e-15).min(1.0 - 1e-15); // Numerical stability
ll += y * p_clamped.ln() + (1.0 - y) * (1.0 - p_clamped).ln();
}
-ll / batch_indices.len() as f64 // Negative log-likelihood
}
}
impl StochasticGradientFunction for LogisticRegressionData {
fn compute_gradient(&mut self, params: &ArrayView1<f64>, batch_indices: &[f64]) -> Array1<f64> {
let mut grad = Array1::zeros(params.len());
for &idx in batch_indices {
let i = (idx as usize) % self.n_samples;
let features = self.features.row(i);
let y = self.labels[i];
let p = self.predict(params, &features);
let error = p - y;
for j in 0..params.len() {
grad[j] += error * features[j];
}
}
grad / batch_indices.len() as f64
}
fn compute_value(&mut self, params: &ArrayView1<f64>, batch_indices: &[f64]) -> f64 {
self.log_likelihood(params, batch_indices)
}
}
fn generate_classification_data(n_samples: usize, n_features: usize) -> (Array2<f64>, Array1<f64>) {
let mut rng = rng();
// Generate features
let mut features = Array2::zeros((n_samples, n_features + 1)); // +1 for bias term
for i in 0..n_samples {
features[[i, 0]] = 1.0; // Bias term
for j in 1..n_features + 1 {
features[[i, j]] = rng.gen_range(-2.0..2.0);
}
}
// Generate true parameters
let true_params: Array1<f64> = (0..n_features + 1)
.map(|_| rng.gen_range(-1.0..1.0))
.collect();
// Generate labels
let labels: Array1<f64> = (0..n_samples)
.map(|i| {
let logit = features.row(i).dot(&true_params);
let p = 1.0 / (1.0 + (-logit).exp());
if rng.gen::<f64>() < p { 1.0 } else { 0.0 }
})
.collect();
(features, labels)
}
fn main() -> Result<(), Box<dyn std::error::Error>> {
// Generate synthetic classification data
let (features, labels) = generate_classification_data(1000, 5);
println!("Generated {} samples with {} features", features.nrows(), features.ncols());
// Split into train/test
let train_size = 800;
let train_features = features.slice(s![0..train_size, ..]).to_owned();
let train_labels = labels.slice(s![0..train_size]).to_owned();
let test_features = features.slice(s![train_size.., ..]).to_owned();
let test_labels = labels.slice(s![train_size..]).to_owned();
// Create logistic regression model
let mut logreg = LogisticRegressionData::new(train_features, train_labels);
// Initial parameters (random initialization)
let mut rng = rng();
let x0: Array1<f64> = (0..features.ncols())
.map(|_| rng.gen_range(-0.1..0.1))
.collect();
// Create data provider
let data_provider = Box::new(InMemoryDataProvider::new(
(0..train_size).map(|i| i as f64).collect()
));
// Training with different optimizers
let optimizers = vec![
("SGD", OptimizerType::SGD),
("Adam", OptimizerType::Adam),
("AdamW", OptimizerType::AdamW),
];
for (name, opt_type) in optimizers {
println!("\nTraining with {}:", name);
let mut model = LogisticRegressionData::new(
logreg.features.clone(),
logreg.labels.clone()
);
let data_provider = Box::new(InMemoryDataProvider::new(
(0..train_size).map(|i| i as f64).collect()
));
let result = match opt_type {
OptimizerType::SGD => {
let options = SGDOptions {
learning_rate: 0.01,
max_iter: 1000,
batch_size: Some(32),
tol: 1e-6,
..Default::default()
};
minimize_sgd(model, x0.clone(), data_provider, options)?
}
OptimizerType::Adam => {
let options = AdamOptions {
learning_rate: 0.001,
max_iter: 1000,
batch_size: Some(32),
tol: 1e-6,
..Default::default()
};
minimize_adam(model, x0.clone(), data_provider, options)?
}
OptimizerType::AdamW => {
let options = AdamWOptions {
learning_rate: 0.001,
max_iter: 1000,
batch_size: Some(32),
weight_decay: 0.01,
tol: 1e-6,
..Default::default()
};
minimize_adamw(model, x0.clone(), data_provider, options)?
}
};
// Evaluate on test set
let mut correct = 0;
for i in 0..test_features.nrows() {
let features = test_features.row(i);
let true_label = test_labels[i];
let z = features.dot(&result.x);
let predicted_prob = 1.0 / (1.0 + (-z).exp());
let predicted_label = if predicted_prob > 0.5 { 1.0 } else { 0.0 };
if predicted_label == true_label {
correct += 1;
}
}
let accuracy = correct as f64 / test_features.nrows() as f64;
println!(" Final loss: {:.6}", result.fun);
println!(" Iterations: {}", result.nit);
println!(" Test accuracy: {:.2}%", accuracy * 100.0);
}
Ok(())
}
#[derive(Clone)]
enum OptimizerType {
SGD,
Adam,
AdamW,
}
```
### 2. Neural Network Training
Train a simple neural network for regression.
```rust
use scirs2_optimize::stochastic::*;
use ndarray::{Array1, Array2, ArrayView1, Axis};
use rand::prelude::*;
struct NeuralNetwork {
features: Array2<f64>,
targets: Array1<f64>,
architecture: Vec<usize>, // [input_size, hidden1, hidden2, ..., output_size]
n_samples: usize,
}
impl NeuralNetwork {
fn new(features: Array2<f64>, targets: Array1<f64>, hidden_sizes: Vec<usize>) -> Self {
let n_samples = features.nrows();
let input_size = features.ncols();
let output_size = 1; // Regression
let mut architecture = vec![input_size];
architecture.extend(hidden_sizes);
architecture.push(output_size);
Self { features, targets, architecture, n_samples }
}
fn get_weight_ranges(&self) -> Vec<(usize, usize)> {
let mut ranges = Vec::new();
let mut start = 0;
for i in 0..self.architecture.len() - 1 {
let weight_count = self.architecture[i] * self.architecture[i + 1];
let bias_count = self.architecture[i + 1];
ranges.push((start, start + weight_count));
ranges.push((start + weight_count, start + weight_count + bias_count));
start += weight_count + bias_count;
}
ranges
}
fn forward(&self, params: &ArrayView1<f64>, input: &ArrayView1<f64>) -> f64 {
let mut activation = input.to_owned();
let ranges = self.get_weight_ranges();
for layer in 0..self.architecture.len() - 1 {
let weight_range = &ranges[layer * 2];
let bias_range = &ranges[layer * 2 + 1];
let weights = params.slice(s![weight_range.0..weight_range.1]);
let biases = params.slice(s![bias_range.0..bias_range.1]);
// Reshape weights to matrix
let weight_matrix = Array2::from_shape_vec(
(self.architecture[layer], self.architecture[layer + 1]),
weights.to_vec(),
).unwrap();
// Linear transformation
let mut new_activation = Array1::zeros(self.architecture[layer + 1]);
for j in 0..self.architecture[layer + 1] {
for i in 0..self.architecture[layer] {
new_activation[j] += activation[i] * weight_matrix[[i, j]];
}
new_activation[j] += biases[j];
}
// Apply activation function (tanh for hidden layers, linear for output)
if layer < self.architecture.len() - 2 {
new_activation.mapv_inplace(|x| x.tanh());
}
activation = new_activation;
}
activation[0] // Single output for regression
}
fn total_parameters(&self) -> usize {
let mut total = 0;
for i in 0..self.architecture.len() - 1 {
total += self.architecture[i] * self.architecture[i + 1]; // Weights
total += self.architecture[i + 1]; // Biases
}
total
}
}
impl StochasticGradientFunction for NeuralNetwork {
fn compute_gradient(&mut self, params: &ArrayView1<f64>, batch_indices: &[f64]) -> Array1<f64> {
let mut grad = Array1::zeros(params.len());
let h = 1e-7; // Finite difference step
// Numerical gradient (for simplicity)
for i in 0..params.len() {
let mut params_plus = params.to_owned();
let mut params_minus = params.to_owned();
params_plus[i] += h;
params_minus[i] -= h;
let loss_plus = self.compute_value(¶ms_plus.view(), batch_indices);
let loss_minus = self.compute_value(¶ms_minus.view(), batch_indices);
grad[i] = (loss_plus - loss_minus) / (2.0 * h);
}
grad
}
fn compute_value(&mut self, params: &ArrayView1<f64>, batch_indices: &[f64]) -> f64 {
let mut loss = 0.0;
for &idx in batch_indices {
let i = (idx as usize) % self.n_samples;
let features = self.features.row(i);
let target = self.targets[i];
let prediction = self.forward(params, &features);
loss += 0.5 * (prediction - target).powi(2); // MSE
}
loss / batch_indices.len() as f64
}
}
// Generate synthetic regression data
fn generate_regression_data(n_samples: usize, noise_level: f64) -> (Array2<f64>, Array1<f64>) {
let mut rng = rng();
// Features: x1, x2
let features = Array2::from_shape_fn((n_samples, 2), |(i, j)| {
rng.gen_range(-2.0..2.0)
});
// Target: y = sin(x1) + cos(x2) + noise
let targets: Array1<f64> = features.outer_iter()
.map(|row| {
let x1 = row[0];
let x2 = row[1];
x1.sin() + x2.cos() + rng.gen_range(-noise_level..noise_level)
})
.collect();
(features, targets)
}
fn main() -> Result<(), Box<dyn std::error::Error>> {
// Generate synthetic data
let (features, targets) = generate_regression_data(1000, 0.1);
println!("Generated {} samples for regression", features.nrows());
// Create neural network: 2 -> 10 -> 5 -> 1
let hidden_sizes = vec![10, 5];
let mut nn = NeuralNetwork::new(features, targets, hidden_sizes);
println!("Neural network architecture: {:?}", nn.architecture);
println!("Total parameters: {}", nn.total_parameters());
// Random initialization
let mut rng = rng();
let x0: Array1<f64> = (0..nn.total_parameters())
.map(|_| rng.gen_range(-0.1..0.1))
.collect();
// Create data provider
let data_provider = Box::new(InMemoryDataProvider::new(
(0..nn.n_samples).map(|i| i as f64).collect()
));
// Train with Adam
let options = AdamOptions {
learning_rate: 0.001,
max_iter: 2000,
batch_size: Some(32),
tol: 1e-6,
..Default::default()
};
let start_time = std::time::Instant::now();
let result = minimize_adam(nn, x0, data_provider, options)?;
let training_time = start_time.elapsed();
println!("\nTraining Results:");
println!(" Final loss: {:.6}", result.fun);
println!(" Iterations: {}", result.nit);
println!(" Training time: {:.2?}", training_time);
println!(" Success: {}", result.success);
// Test the trained network
let test_features = Array2::from_shape_vec((5, 2), vec![
0.0, 0.0, // sin(0) + cos(0) = 1
1.0, 0.0, // sin(1) + cos(0) ≈ 1.84
0.0, 1.0, // sin(0) + cos(1) ≈ 0.54
1.5, 1.5, // sin(1.5) + cos(1.5) ≈ 1.07
-1.0, -1.0, // sin(-1) + cos(-1) ≈ -0.30
])?;
println!("\nTest Predictions:");
for (i, row) in test_features.outer_iter().enumerate() {
let x1 = row[0];
let x2 = row[1];
let true_value = x1.sin() + x2.cos();
// Create a dummy NN for prediction (we need the architecture)
let dummy_features = Array2::zeros((1, 2));
let dummy_targets = Array1::zeros(1);
let dummy_nn = NeuralNetwork::new(dummy_features, dummy_targets, vec![10, 5]);
let predicted = dummy_nn.forward(&result.x.view(), &row);
println!(" x=[{:.1}, {:.1}]: true={:.3}, predicted={:.3}, error={:.3}",
x1, x2, true_value, predicted, (true_value - predicted).abs());
}
Ok(())
}
```
## Engineering Applications
### 1. Parameter Identification in Dynamical Systems
Identify parameters of a damped oscillator from experimental data.
```rust
use scirs2_optimize::least_squares::*;
use ndarray::{Array1, Array2};
use std::f64::consts::PI;
// Model: x(t) = A * exp(-ζωt) * cos(ωd*t + φ)
// where ωd = ω * sqrt(1 - ζ²) for underdamped case
struct DampedOscillator {
time_data: Vec<f64>,
position_data: Vec<f64>,
}
impl DampedOscillator {
fn new(time_data: Vec<f64>, position_data: Vec<f64>) -> Self {
assert_eq!(time_data.len(), position_data.len());
Self { time_data, position_data }
}
// Parameters: [A, ω, ζ, φ]
fn model(&self, params: &[f64], t: f64) -> f64 {
let (a, omega, zeta, phi) = (params[0], params[1], params[2], params[3]);
if zeta >= 1.0 {
// Overdamped case
let r1 = -omega * (zeta + (zeta.powi(2) - 1.0).sqrt());
let r2 = -omega * (zeta - (zeta.powi(2) - 1.0).sqrt());
a * ((r1 * t).exp() + (r2 * t).exp()) // Simplified for demonstration
} else {
// Underdamped case
let omega_d = omega * (1.0 - zeta.powi(2)).sqrt();
a * (-zeta * omega * t).exp() * (omega_d * t + phi).cos()
}
}
fn residuals(&self, params: &[f64]) -> Array1<f64> {
let mut residuals = Array1::zeros(self.time_data.len());
for (i, (&t, &x_measured)) in self.time_data.iter()
.zip(self.position_data.iter())
.enumerate() {
let x_predicted = self.model(params, t);
residuals[i] = x_measured - x_predicted;
}
residuals
}
}
// Generate synthetic experimental data
fn generate_oscillator_data(
true_params: &[f64],
t_max: f64,
n_points: usize,
noise_level: f64
) -> (Vec<f64>, Vec<f64>) {
use rand::prelude::*;
let mut rng = rng();
let dt = t_max / (n_points - 1) as f64;
let time_data: Vec<f64> = (0..n_points).map(|i| i as f64 * dt).collect();
let oscillator = DampedOscillator::new(vec![], vec![]);
let position_data: Vec<f64> = time_data.iter()
.map(|&t| {
let true_value = oscillator.model(true_params, t);
true_value + rng.gen_range(-noise_level..noise_level)
})
.collect();
(time_data, position_data)
}
fn main() -> Result<(), Box<dyn std::error::Error>> {
// True parameters: A=2.0, ω=5.0, ζ=0.1, φ=π/4
let true_params = vec![2.0, 5.0, 0.1, PI / 4.0];
// Generate synthetic experimental data
let (time_data, position_data) = generate_oscillator_data(&true_params, 2.0, 100, 0.05);
println!("Generated {} data points for parameter identification", time_data.len());
println!("True parameters: A={:.2}, ω={:.2}, ζ={:.3}, φ={:.3}",
true_params[0], true_params[1], true_params[2], true_params[3]);
let oscillator = DampedOscillator::new(time_data, position_data);
// Residual function for least squares
let residual_fn = |params: &[f64], _data: &[f64]| -> Array1<f64> {
oscillator.residuals(params)
};
// Initial guess (perturbed true values)
let x0 = Array1::from_vec(vec![1.5, 4.0, 0.2, 0.5]);
let dummy_data = Array1::zeros(0);
println!("\nInitial guess: A={:.2}, ω={:.2}, ζ={:.3}, φ={:.3}",
x0[0], x0[1], x0[2], x0[3]);
// Parameter identification using different methods
let methods = vec![
("Levenberg-Marquardt", Method::LevenbergMarquardt),
("Trust Region Reflective", Method::TrustRegionReflective),
];
for (name, method) in methods {
println!("\n{} Results:", name);
let result = least_squares(
residual_fn,
&x0,
method,
None, // Auto-compute Jacobian
&dummy_data,
None
)?;
println!(" Identified parameters:");
println!(" A = {:.4} (true: {:.4}, error: {:.2}%)",
result.x[0], true_params[0],
100.0 * (result.x[0] - true_params[0]).abs() / true_params[0]);
println!(" ω = {:.4} (true: {:.4}, error: {:.2}%)",
result.x[1], true_params[1],
100.0 * (result.x[1] - true_params[1]).abs() / true_params[1]);
println!(" ζ = {:.4} (true: {:.4}, error: {:.2}%)",
result.x[2], true_params[2],
100.0 * (result.x[2] - true_params[2]).abs() / true_params[2]);
println!(" φ = {:.4} (true: {:.4}, error: {:.2}%)",
result.x[3], true_params[3],
100.0 * (result.x[3] - true_params[3]).abs() / true_params[3]);
println!(" Final residual norm: {:.2e}", result.fun);
println!(" Iterations: {}", result.nit);
println!(" Success: {}", result.success);
// Compute R² coefficient of determination
let residuals = oscillator.residuals(&result.x.to_vec());
let ss_res: f64 = residuals.iter().map(|r| r.powi(2)).sum();
let y_mean: f64 = oscillator.position_data.iter().sum::<f64>() / oscillator.position_data.len() as f64;
let ss_tot: f64 = oscillator.position_data.iter()
.map(|y| (y - y_mean).powi(2))
.sum();
let r_squared = 1.0 - ss_res / ss_tot;
println!(" R² = {:.6}", r_squared);
}
Ok(())
}
```
### 2. Optimal Control Problem
Solve a trajectory optimization problem for a robotic arm.
```rust
use scirs2_optimize::prelude::*;
use ndarray::{Array1, Array2};
use std::f64::consts::PI;
struct RobotArmOptimization {
n_timesteps: usize,
dt: f64,
target_position: [f64; 2],
arm_lengths: [f64; 2],
q_weight: f64, // State cost weight
r_weight: f64, // Control cost weight
qf_weight: f64, // Terminal cost weight
}
impl RobotArmOptimization {
fn new(n_timesteps: usize, dt: f64) -> Self {
Self {
n_timesteps,
dt,
target_position: [1.0, 1.0],
arm_lengths: [1.0, 1.0],
q_weight: 1.0,
r_weight: 0.1,
qf_weight: 100.0,
}
}
// Forward kinematics: compute end-effector position
fn forward_kinematics(&self, angles: &[f64]) -> [f64; 2] {
let (q1, q2) = (angles[0], angles[1]);
let (l1, l2) = (self.arm_lengths[0], self.arm_lengths[1]);
let x = l1 * q1.cos() + l2 * (q1 + q2).cos();
let y = l1 * q1.sin() + l2 * (q1 + q2).sin();
[x, y]
}
// Dynamics: simple double integrator model
fn dynamics(&self, state: &[f64], control: &[f64]) -> [f64; 4] {
// State: [q1, q2, q1_dot, q2_dot]
// Control: [tau1, tau2] (torques)
let q1_dot = state[2];
let q2_dot = state[3];
let q1_ddot = control[0]; // Simplified dynamics
let q2_ddot = control[1];
[q1_dot, q2_dot, q1_ddot, q2_ddot]
}
// Integrate dynamics using Euler method
fn integrate_step(&self, state: &[f64], control: &[f64]) -> [f64; 4] {
let state_dot = self.dynamics(state, control);
[
state[0] + self.dt * state_dot[0],
state[1] + self.dt * state_dot[1],
state[2] + self.dt * state_dot[2],
state[3] + self.dt * state_dot[3],
]
}
// Cost function for trajectory optimization
fn trajectory_cost(&self, decision_vars: &ArrayView1<f64>) -> f64 {
// Decision variables: [u0, u1, ..., u_{N-1}] where each u_i = [tau1_i, tau2_i]
// Initial state is fixed: [0, 0, 0, 0]
let mut state = [0.0, 0.0, 0.0, 0.0]; // [q1, q2, q1_dot, q2_dot]
let mut total_cost = 0.0;
// Simulate trajectory and accumulate cost
for k in 0..self.n_timesteps {
let control = [decision_vars[2 * k], decision_vars[2 * k + 1]];
// Running cost: state deviation + control effort
let end_effector = self.forward_kinematics(&state[0..2]);
let position_error = [
end_effector[0] - self.target_position[0],
end_effector[1] - self.target_position[1],
];
let state_cost = self.q_weight * (
position_error[0].powi(2) + position_error[1].powi(2) +
0.1 * (state[2].powi(2) + state[3].powi(2)) // Velocity penalty
);
let control_cost = self.r_weight * (control[0].powi(2) + control[1].powi(2));
total_cost += self.dt * (state_cost + control_cost);
// Integrate dynamics
state = self.integrate_step(&state, &control);
}
// Terminal cost
let final_end_effector = self.forward_kinematics(&state[0..2]);
let final_error = [
final_end_effector[0] - self.target_position[0],
final_end_effector[1] - self.target_position[1],
];
let terminal_cost = self.qf_weight * (
final_error[0].powi(2) + final_error[1].powi(2)
);
total_cost + terminal_cost
}
}
fn main() -> Result<(), Box<dyn std::error::Error>> {
let n_timesteps = 50;
let dt = 0.02; // 20ms timesteps
let robot = RobotArmOptimization::new(n_timesteps, dt);
println!("Robot Arm Trajectory Optimization");
println!(" Target position: [{:.2}, {:.2}]",
robot.target_position[0], robot.target_position[1]);
println!(" Time horizon: {:.1}s ({} timesteps)", n_timesteps as f64 * dt, n_timesteps);
// Decision variables: torques for each timestep [tau1_0, tau2_0, tau1_1, tau2_1, ...]
let n_vars = 2 * n_timesteps;
// Initial guess: zero torques
let x0 = Array1::zeros(n_vars);
let objective = |x: &ArrayView1<f64>| robot.trajectory_cost(x);
// Add bounds on torques: -10 to 10 Nm
let bounds = Bounds::new(&vec![(Some(-10.0), Some(10.0)); n_vars]);
let mut options = Options::default();
options.bounds = Some(bounds);
options.max_iter = 1000;
options.gtol = 1e-6;
println!("\nOptimizing trajectory...");
let start_time = std::time::Instant::now();
let result = minimize(objective, &x0, Method::LBFGS, Some(options))?;
let optimization_time = start_time.elapsed();
println!("\nOptimization Results:");
println!(" Final cost: {:.6}", result.fun);
println!(" Iterations: {}", result.nit);
println!(" Function evaluations: {}", result.nfev);
println!(" Optimization time: {:.2?}", optimization_time);
println!(" Success: {}", result.success);
// Simulate the optimal trajectory
println!("\nSimulating optimal trajectory:");
let mut state = [0.0, 0.0, 0.0, 0.0];
println!("Time q1 q2 x_ee y_ee Distance to target");
for k in 0..n_timesteps.min(10) { // Show first 10 timesteps
let control = [result.x[2 * k], result.x[2 * k + 1]];
let end_effector = robot.forward_kinematics(&state[0..2]);
let distance = ((end_effector[0] - robot.target_position[0]).powi(2) +
(end_effector[1] - robot.target_position[1]).powi(2)).sqrt();
println!("{:.2} {:.3} {:.3} {:.3} {:.3} {:.3}",
k as f64 * dt, state[0], state[1],
end_effector[0], end_effector[1], distance);
state = robot.integrate_step(&state, &control);
}
// Final position
let final_end_effector = robot.forward_kinematics(&state[0..2]);
let final_distance = ((final_end_effector[0] - robot.target_position[0]).powi(2) +
(final_end_effector[1] - robot.target_position[1]).powi(2)).sqrt();
println!("...");
println!("Final end-effector position: [{:.4}, {:.4}]",
final_end_effector[0], final_end_effector[1]);
println!("Final distance to target: {:.4}", final_distance);
// Analyze control effort
let max_torque = result.x.iter().map(|&u| u.abs()).fold(0.0, f64::max);
let avg_torque = result.x.iter().map(|&u| u.abs()).sum::<f64>() / result.x.len() as f64;
println!("\nControl Analysis:");
println!(" Maximum torque: {:.2} Nm", max_torque);
println!(" Average torque magnitude: {:.2} Nm", avg_torque);
Ok(())
}
```
## Financial Applications
### 1. Portfolio Optimization with Risk Management
Modern portfolio theory with additional constraints.
```rust
use scirs2_optimize::{constrained::*, prelude::*};
use ndarray::{Array1, Array2};
struct PortfolioOptimizer {
expected_returns: Array1<f64>,
covariance_matrix: Array2<f64>,
n_assets: usize,
}
impl PortfolioOptimizer {
fn new(expected_returns: Array1<f64>, covariance_matrix: Array2<f64>) -> Self {
let n_assets = expected_returns.len();
assert_eq!(covariance_matrix.dim(), (n_assets, n_assets));
Self { expected_returns, covariance_matrix, n_assets }
}
// Mean-variance utility: μ'w - (λ/2)w'Σw
fn utility(&self, weights: &ArrayView1<f64>, risk_aversion: f64) -> f64 {
let expected_return = self.expected_returns.dot(weights);
let portfolio_variance = weights.dot(&self.covariance_matrix.dot(weights));
expected_return - 0.5 * risk_aversion * portfolio_variance
}
// Portfolio variance
fn portfolio_variance(&self, weights: &ArrayView1<f64>) -> f64 {
weights.dot(&self.covariance_matrix.dot(weights))
}
// Portfolio return
fn portfolio_return(&self, weights: &ArrayView1<f64>) -> f64 {
self.expected_returns.dot(weights)
}
// Maximum drawdown constraint (simplified)
fn max_drawdown_constraint(&self, weights: &ArrayView1<f64>, max_dd: f64) -> f64 {
// Simplified: assume drawdown is proportional to volatility
let volatility = self.portfolio_variance(weights).sqrt();
max_dd - 2.0 * volatility // Rough approximation
}
// Sector exposure constraint
fn sector_constraint(&self, weights: &ArrayView1<f64>, sector_indices: &[usize], max_exposure: f64) -> f64 {
let sector_weight: f64 = sector_indices.iter().map(|&i| weights[i]).sum();
max_exposure - sector_weight
}
}
fn main() -> Result<(), Box<dyn std::error::Error>> {
// Example data: 5 assets with different characteristics
let asset_names = vec!["Stocks", "Bonds", "REITs", "Commodities", "Cash"];
let n_assets = asset_names.len();
// Annual expected returns (%)
let expected_returns = Array1::from_vec(vec![0.10, 0.04, 0.08, 0.06, 0.02]);
// Covariance matrix (annual)
let covariance_data = vec![
0.0400, 0.0012, 0.0080, 0.0100, 0.0000, // Stocks
0.0012, 0.0025, 0.0010, 0.0020, 0.0000, // Bonds
0.0080, 0.0010, 0.0225, 0.0050, 0.0000, // REITs
0.0100, 0.0020, 0.0050, 0.0625, 0.0000, // Commodities
0.0000, 0.0000, 0.0000, 0.0000, 0.0001, // Cash
];
let covariance_matrix = Array2::from_shape_vec((n_assets, n_assets), covariance_data)?;
let portfolio = PortfolioOptimizer::new(expected_returns, covariance_matrix);
println!("Portfolio Optimization with Risk Management");
println!("Assets: {:?}", asset_names);
println!("Expected returns: {:?}", portfolio.expected_returns);
// Risk aversion parameter
let risk_aversion = 5.0;
// Objective: maximize utility = minimize negative utility
let objective = |w: &ArrayView1<f64>| -> f64 {
-portfolio.utility(w, risk_aversion)
};
// Constraints
let mut constraints = Vec::new();
// 1. Weights sum to 1
constraints.push(Constraint {
fun: Box::new(|w: &ArrayView1<f64>| w.sum() - 1.0),
constraint_type: ConstraintType::Equality,
jac: None,
});
// 2. Minimum return constraint (6% annually)
let min_return = 0.06;
constraints.push(Constraint {
fun: Box::new(move |w: &ArrayView1<f64>| {
portfolio.portfolio_return(w) - min_return
}),
constraint_type: ConstraintType::Inequality,
jac: None,
});
// 3. Maximum volatility constraint (15% annually)
let max_volatility = 0.15;
constraints.push(Constraint {
fun: Box::new(move |w: &ArrayView1<f64>| {
max_volatility.powi(2) - portfolio.portfolio_variance(w)
}),
constraint_type: ConstraintType::Inequality,
jac: None,
});
// 4. Sector constraints: max 60% in risky assets (stocks + REITs + commodities)
let risky_indices = vec![0, 2, 3]; // Stocks, REITs, Commodities
let max_risky_exposure = 0.60;
constraints.push(Constraint {
fun: Box::new(move |w: &ArrayView1<f64>| {
portfolio.sector_constraint(w, &risky_indices, max_risky_exposure)
}),
constraint_type: ConstraintType::Inequality,
jac: None,
});
// Bounds: 0 <= w_i <= 0.4 (max 40% in any single asset)
let bounds = Bounds::new(&vec![(Some(0.0), Some(0.4)); n_assets]);
// Initial guess: equal weights
let x0 = Array1::from_elem(n_assets, 1.0 / n_assets as f64);
let mut options = Options::default();
options.bounds = Some(bounds);
options.max_iter = 1000;
println!("\nOptimizing portfolio...");
let start_time = std::time::Instant::now();
let result = minimize_constrained(objective, &x0, &constraints, Method::SLSQP, Some(options))?;
let optimization_time = start_time.elapsed();
println!("\nOptimization Results:");
println!(" Success: {}", result.success);
println!(" Iterations: {}", result.nit);
println!(" Optimization time: {:.2?}", optimization_time);
// Analyze optimal portfolio
println!("\nOptimal Portfolio Allocation:");
for (i, &weight) in result.x.iter().enumerate() {
println!(" {}: {:.1}%", asset_names[i], weight * 100.0);
}
let portfolio_return = portfolio.portfolio_return(&result.x.view());
let portfolio_vol = portfolio.portfolio_variance(&result.x.view()).sqrt();
let sharpe_ratio = (portfolio_return - 0.02) / portfolio_vol; // Assume 2% risk-free rate
println!("\nPortfolio Characteristics:");
println!(" Expected return: {:.2}%", portfolio_return * 100.0);
println!(" Volatility: {:.2}%", portfolio_vol * 100.0);
println!(" Sharpe ratio: {:.3}", sharpe_ratio);
// Verify constraints
println!("\nConstraint Verification:");
println!(" Sum of weights: {:.6}", result.x.sum());
println!(" Return constraint: {:.2}% >= {:.2}%",
portfolio_return * 100.0, min_return * 100.0);
println!(" Volatility constraint: {:.2}% <= {:.2}%",
portfolio_vol * 100.0, max_volatility * 100.0);
let risky_allocation: f64 = risky_indices.iter().map(|&i| result.x[i]).sum();
println!(" Risky assets allocation: {:.1}% <= {:.1}%",
risky_allocation * 100.0, max_risky_exposure * 100.0);
// Risk analysis
println!("\nRisk Analysis:");
// Value at Risk (simplified normal approximation)
let confidence_levels = vec![0.95, 0.99];
for &confidence in &confidence_levels {
let z_score = match confidence {
0.95 => 1.645,
0.99 => 2.326,
_ => 1.960,
};
let var = z_score * portfolio_vol;
println!(" VaR ({}%): {:.2}%", (confidence * 100.0) as i32, var * 100.0);
}
// Component contributions to risk
println!("\nRisk Contributions:");
let total_risk = portfolio_vol.powi(2);
for (i, &weight) in result.x.iter().enumerate() {
let marginal_risk = portfolio.covariance_matrix.row(i).dot(&result.x);
let contribution = weight * marginal_risk / total_risk;
println!(" {}: {:.1}%", asset_names[i], contribution * 100.0);
}
Ok(())
}
```
## Scientific Computing
### 1. Protein Folding Energy Minimization
Simplified protein folding using optimization.
```rust
use scirs2_optimize::{global::*, prelude::*};
use ndarray::Array1;
use std::f64::consts::PI;
struct ProteinFolder {
n_residues: usize,
bond_length: f64,
amino_acid_types: Vec<u8>, // 0-19 for 20 amino acids
}
impl ProteinFolder {
fn new(sequence: &str) -> Self {
let amino_acids = "ACDEFGHIKLMNPQRSTVWY";
let amino_acid_types: Vec<u8> = sequence.chars()
.map(|c| amino_acids.find(c).unwrap_or(0) as u8)
.collect();
Self {
n_residues: sequence.len(),
bond_length: 3.8, // Angstroms
amino_acid_types,
}
}
// Convert angles to 3D coordinates
fn angles_to_coordinates(&self, angles: &ArrayView1<f64>) -> Vec<[f64; 3]> {
let mut coords = vec![[0.0; 3]; self.n_residues];
// First residue at origin
coords[0] = [0.0, 0.0, 0.0];
if self.n_residues > 1 {
// Second residue along x-axis
coords[1] = [self.bond_length, 0.0, 0.0];
}
// Build the rest using spherical coordinates
for i in 2..self.n_residues {
let phi = angles[2 * (i - 2)]; // Dihedral angle
let theta = angles[2 * (i - 2) + 1]; // Bond angle
// Previous two residues
let p1 = coords[i - 1];
let p2 = coords[i - 2];
// Direction vector from p2 to p1
let v = [p1[0] - p2[0], p1[1] - p2[1], p1[2] - p2[2]];
let v_len = (v[0].powi(2) + v[1].powi(2) + v[2].powi(2)).sqrt();
let v_norm = [v[0] / v_len, v[1] / v_len, v[2] / v_len];
// New position using spherical coordinates
let x = p1[0] + self.bond_length * (theta.cos() * v_norm[0] -
theta.sin() * phi.cos() * v_norm[1] -
theta.sin() * phi.sin() * v_norm[2]);
let y = p1[1] + self.bond_length * (theta.cos() * v_norm[1] +
theta.sin() * phi.cos() * v_norm[0]);
let z = p1[2] + self.bond_length * (theta.cos() * v_norm[2] +
theta.sin() * phi.sin() * v_norm[0]);
coords[i] = [x, y, z];
}
coords
}
// Simplified energy function
fn energy(&self, angles: &ArrayView1<f64>) -> f64 {
let coords = self.angles_to_coordinates(angles);
let mut total_energy = 0.0;
// 1. Bond length energy (harmonic potential)
for i in 1..self.n_residues {
let dist = self.distance(&coords[i - 1], &coords[i]);
let bond_energy = 100.0 * (dist - self.bond_length).powi(2);
total_energy += bond_energy;
}
// 2. Lennard-Jones potential for non-bonded interactions
for i in 0..self.n_residues {
for j in i + 2..self.n_residues { // Skip adjacent residues
let dist = self.distance(&coords[i], &coords[j]);
let sigma = self.get_interaction_distance(i, j);
let epsilon = self.get_interaction_strength(i, j);
if dist > 0.1 { // Avoid division by zero
let r6 = (sigma / dist).powi(6);
let lj_energy = 4.0 * epsilon * (r6.powi(2) - r6);
total_energy += lj_energy;
}
}
}
// 3. Electrostatic interactions (simplified)
for i in 0..self.n_residues {
for j in i + 1..self.n_residues {
let dist = self.distance(&coords[i], &coords[j]);
let charge_i = self.get_charge(i);
let charge_j = self.get_charge(j);
if dist > 0.1 {
let electrostatic = 332.0 * charge_i * charge_j / dist; // kcal/mol
total_energy += electrostatic;
}
}
}
// 4. Secondary structure preferences
for i in 0..self.n_residues.saturating_sub(3) {
let alpha_helix_penalty = self.alpha_helix_penalty(&angles, i);
let beta_sheet_penalty = self.beta_sheet_penalty(&coords, i);
total_energy += alpha_helix_penalty + beta_sheet_penalty;
}
total_energy
}
fn distance(&self, p1: &[f64; 3], p2: &[f64; 3]) -> f64 {
((p1[0] - p2[0]).powi(2) + (p1[1] - p2[1]).powi(2) + (p1[2] - p2[2]).powi(2)).sqrt()
}
fn get_interaction_distance(&self, i: usize, j: usize) -> f64 {
// Simplified: based on amino acid types
let type_i = self.amino_acid_types[i];
let type_j = self.amino_acid_types[j];
3.5 + 0.1 * (type_i as f64 + type_j as f64) / 40.0 // 3.5-4.0 Angstroms
}
fn get_interaction_strength(&self, i: usize, j: usize) -> f64 {
// Simplified: hydrophobic interactions
let hydrophobic = [0, 5, 8, 11, 15, 17, 18, 19]; // Indices of hydrophobic amino acids
let type_i = self.amino_acid_types[i];
let type_j = self.amino_acid_types[j];
let is_hydrophobic_i = hydrophobic.contains(&type_i);
let is_hydrophobic_j = hydrophobic.contains(&type_j);
match (is_hydrophobic_i, is_hydrophobic_j) {
(true, true) => 1.0, // Strong hydrophobic interaction
(true, false) | (false, true) => 0.5, // Weak interaction
(false, false) => 0.2, // Hydrophilic interaction
}
}
fn get_charge(&self, i: usize) -> f64 {
// Simplified charges for amino acids
match self.amino_acid_types[i] {
3 => -1.0, // D (Asp)
4 => -1.0, // E (Glu)
9 => 1.0, // K (Lys)
14 => 1.0, // R (Arg)
6 => 0.5, // H (His)
_ => 0.0, // Neutral
}
}
fn alpha_helix_penalty(&self, angles: &ArrayView1<f64>, start_idx: usize) -> f64 {
// Ideal alpha helix: phi ≈ -60°, psi ≈ -45°
let ideal_phi = -60.0 * PI / 180.0;
let ideal_psi = -45.0 * PI / 180.0;
let mut penalty = 0.0;
for i in 0..4.min(self.n_residues - start_idx - 2) {
let idx = start_idx + i;
if 2 * idx + 1 < angles.len() {
let phi = angles[2 * idx];
let psi = angles[2 * idx + 1];
let phi_diff = (phi - ideal_phi).abs();
let psi_diff = (psi - ideal_psi).abs();
penalty += 0.1 * (phi_diff.powi(2) + psi_diff.powi(2));
}
}
penalty
}
fn beta_sheet_penalty(&self, coords: &[[f64; 3]], start_idx: usize) -> f64 {
// Simplified: penalize extended conformations
if start_idx + 2 >= coords.len() {
return 0.0;
}
let d1 = self.distance(&coords[start_idx], &coords[start_idx + 1]);
let d2 = self.distance(&coords[start_idx + 1], &coords[start_idx + 2]);
let d3 = self.distance(&coords[start_idx], &coords[start_idx + 2]);
// Beta sheet preference: extended conformation
let extension = d3 / (d1 + d2);
if extension > 0.8 {
-0.5 // Reward extended conformations
} else {
0.0
}
}
}
fn main() -> Result<(), Box<dyn std::error::Error>> {
// Example protein sequence (simplified)
let sequence = "MKKLKKLLKKLLKK"; // Alternating hydrophobic/hydrophilic
let folder = ProteinFolder::new(sequence);
println!("Protein Folding Optimization");
println!("Sequence: {}", sequence);
println!("Length: {} residues", folder.n_residues);
// Two angles per residue (except first two): phi and psi
let n_angles = 2 * (folder.n_residues - 2);
println!("Degrees of freedom: {} angles", n_angles);
// Energy function
let energy_func = |angles: &ArrayView1<f64>| folder.energy(angles);
// Bounds: angles in [-π, π]
let bounds: Vec<(f64, f64)> = (0..n_angles).map(|_| (-PI, PI)).collect();
println!("\nOptimizing protein conformation...");
// Try different global optimization methods
let methods = vec![
("Differential Evolution", "DE"),
("Simulated Annealing", "SA"),
("Basin-hopping", "BH"),
];
let mut best_energy = f64::INFINITY;
let mut best_conformation = Array1::zeros(n_angles);
for (method_name, method_code) in methods {
println!("\n{} optimization:", method_name);
let start_time = std::time::Instant::now();
let result = match method_code {
"DE" => {
let options = DifferentialEvolutionOptions {
population_size: 50,
max_iter: 1000,
f: 0.8,
cr: 0.9,
tol: 1e-6,
..Default::default()
};
differential_evolution(energy_func, &bounds, Some(options))?
}
"SA" => {
let options = SimulatedAnnealingOptions {
max_iter: 10000,
initial_temp: 1000.0,
cooling_rate: 0.95,
min_temp: 1e-3,
..Default::default()
};
simulated_annealing(energy_func, &bounds, Some(options))?
}
"BH" => {
// Start from random point
use rand::prelude::*;
let mut rng = rng();
let x0: Array1<f64> = (0..n_angles)
.map(|_| rng.gen_range(-PI..PI))
.collect();
let options = BasinHoppingOptions {
n_iter: 100,
temperature: 10.0,
step_size: 0.5,
..Default::default()
};
basinhopping(energy_func, &x0, Some(options))?
}
_ => unreachable!(),
};
let optimization_time = start_time.elapsed();
println!(" Final energy: {:.2} kcal/mol", result.fun);
println!(" Iterations: {}", result.nit);
println!(" Time: {:.2?}", optimization_time);
println!(" Success: {}", result.success);
if result.fun < best_energy {
best_energy = result.fun;
best_conformation = result.x.clone();
println!(" *** New best conformation! ***");
}
}
println!("\nBest conformation analysis:");
println!(" Final energy: {:.2} kcal/mol", best_energy);
// Analyze the best structure
let coords = folder.angles_to_coordinates(&best_conformation.view());
// Radius of gyration
let center = coords.iter().fold([0.0, 0.0, 0.0], |acc, coord| {
[acc[0] + coord[0], acc[1] + coord[1], acc[2] + coord[2]]
});
let center = [
center[0] / coords.len() as f64,
center[1] / coords.len() as f64,
center[2] / coords.len() as f64,
];
let rg_squared = coords.iter()
.map(|coord| {
(coord[0] - center[0]).powi(2) +
(coord[1] - center[1]).powi(2) +
(coord[2] - center[2]).powi(2)
})
.sum::<f64>() / coords.len() as f64;
let radius_of_gyration = rg_squared.sqrt();
println!(" Radius of gyration: {:.2} Å", radius_of_gyration);
// End-to-end distance
let end_to_end = folder.distance(&coords[0], &coords[coords.len() - 1]);
println!(" End-to-end distance: {:.2} Å", end_to_end);
// Compactness
let max_extent = coords.iter()
.flat_map(|coord| coord.iter())
.fold(0.0, |acc, &x| acc.max(x.abs()));
println!(" Maximum extent: {:.2} Å", max_extent);
// Secondary structure analysis (simplified)
println!("\nSecondary structure angles (degrees):");
for i in 0..5.min(n_angles / 2) {
let phi = best_conformation[2 * i] * 180.0 / PI;
let psi = best_conformation[2 * i + 1] * 180.0 / PI;
println!(" Residue {}: φ = {:.1}°, ψ = {:.1}°", i + 3, phi, psi);
}
Ok(())
}
```
This comprehensive examples document demonstrates the versatility and power of the `scirs2-optimize` library across various domains. Each example includes complete, runnable code with detailed explanations of the problem setup, optimization approach, and result analysis.