diff --git a/src/constraint.rs b/src/constraint.rs
index 4b8f5ba..d97dc38 100644
--- a/src/constraint.rs
+++ b/src/constraint.rs
@@ -1,10 +1,15 @@
 use crate::ffi;
+use crate::scip::ScipPtr;
+use std::rc::Rc;
 
 /// A constraint in an optimization problem.
 #[derive(Debug)]
+#[allow(dead_code)]
 pub struct Constraint {
     /// A pointer to the underlying `SCIP_CONS` C struct.
     pub(crate) raw: *mut ffi::SCIP_CONS,
+    /// A reference to the SCIP instance that owns this constraint (to prevent freeing the model while the constraint is live).
+    pub(crate) scip: Rc<ScipPtr>,
 }
 
 impl Constraint {
@@ -22,3 +27,24 @@ impl Constraint {
         }
     }
 }
+
+#[cfg(test)]
+mod tests {
+    use crate::prelude::*;
+
+    #[test]
+    fn test_constraint_mem_safety() {
+        // Create model
+        let mut model = Model::new()
+            .hide_output()
+            .include_default_plugins()
+            .create_prob("test")
+            .set_obj_sense(ObjSense::Maximize);
+
+        let x1 = model.add_var(0., f64::INFINITY, 3., "x1", VarType::Integer);
+        let cons = model.add_cons(vec![x1], &[1.], 4., 4., "cons");
+        drop(model);
+
+        assert_eq!(cons.name(), "cons");
+    }
+}
diff --git a/src/model.rs b/src/model.rs
index d0bff46..efc0186 100644
--- a/src/model.rs
+++ b/src/model.rs
@@ -106,7 +106,18 @@ impl Model<PluginsIncluded> {
         let scip = self.scip.clone();
         scip.read_prob(filename)?;
         let vars = Rc::new(RefCell::new(self.scip.vars()));
-        let conss = Rc::new(RefCell::new(self.scip.conss()));
+        let conss = Rc::new(RefCell::new(
+            self.scip
+                .conss()
+                .iter()
+                .map(|c| {
+                    Rc::new(Constraint {
+                        raw: *c,
+                        scip: self.scip.clone(),
+                    })
+                })
+                .collect(),
+        ));
         let new_model = Model {
             scip: self.scip,
             state: ProblemCreated { vars, conss },
@@ -835,7 +846,12 @@ macro_rules! impl_ProblemOrSolving {
                         name,
                     )
                     .expect("Failed to create constraint in state ProblemCreated");
-                let cons = Rc::new(cons);
+                let cons = Rc::new(
+                    Constraint {
+                        raw: cons,
+                        scip: self.scip.clone(),
+                    }
+                );
                 self.state.conss.borrow_mut().push(cons.clone());
                 cons
             }
@@ -872,7 +888,10 @@ macro_rules! impl_ProblemOrSolving {
                     .scip
                     .create_cons(vars, coefs, lhs, rhs, name)
                     .expect("Failed to create constraint in state ProblemCreated");
-                let cons = Rc::new(cons);
+                let cons = Rc::new( Constraint {
+                    raw: cons,
+                    scip: self.scip.clone(),
+                });
                 self.state.conss.borrow_mut().push(cons.clone());
                 cons
             }
@@ -897,7 +916,10 @@ macro_rules! impl_ProblemOrSolving {
                     .scip
                     .create_cons_set_part(vars, name)
                     .expect("Failed to add constraint set partition in state ProblemCreated");
-                let cons = Rc::new(cons);
+                let cons = Rc::new( Constraint {
+                    raw: cons,
+                    scip: self.scip.clone(),
+                });
                 self.state.conss.borrow_mut().push(cons.clone());
                 cons
             }
@@ -922,7 +944,10 @@ macro_rules! impl_ProblemOrSolving {
                     .scip
                     .create_cons_set_cover(vars, name)
                     .expect("Failed to add constraint set cover in state ProblemCreated");
-                let cons = Rc::new(cons);
+                let cons = Rc::new( Constraint {
+                    raw: cons,
+                    scip: self.scip.clone(),
+                });
                 self.state.conss.borrow_mut().push(cons.clone());
                 cons
             }
@@ -947,7 +972,10 @@ macro_rules! impl_ProblemOrSolving {
                     .scip
                     .create_cons_set_pack(vars, name)
                     .expect("Failed to add constraint set packing in state ProblemCreated");
-                let cons = Rc::new(cons);
+                let cons = Rc::new( Constraint {
+                    raw: cons,
+                    scip: self.scip.clone(),
+                });
                 self.state.conss.borrow_mut().push(cons.clone());
                 cons
             }
@@ -972,7 +1000,10 @@ macro_rules! impl_ProblemOrSolving {
                     .scip
                     .create_cons_cardinality(vars, cardinality, name)
                     .expect("Failed to add cardinality constraint");
-                let cons = Rc::new(cons);
+                let cons = Rc::new( Constraint {
+                    raw: cons,
+                    scip: self.scip.clone(),
+                });
                 self.state.conss.borrow_mut().push(cons.clone());
                 cons
             }
@@ -1008,7 +1039,10 @@ macro_rules! impl_ProblemOrSolving {
                     .scip
                     .create_cons_indicator(bin_var, vars, coefs, rhs, name)
                     .expect("Failed to create constraint in state ProblemCreated");
-                let cons = Rc::new(cons);
+                let cons = Rc::new( Constraint {
+                    raw: cons,
+                    scip: self.scip.clone(),
+                });
                 self.state.conss.borrow_mut().push(cons.clone());
                 cons
             }
diff --git a/src/scip.rs b/src/scip.rs
index f515f1b..d13814a 100644
--- a/src/scip.rs
+++ b/src/scip.rs
@@ -6,7 +6,7 @@ use crate::{
 };
 use crate::{scip_call, HeurTiming, Heuristic};
 use core::panic;
-use scip_sys::SCIP_SOL;
+use scip_sys::{SCIP_Cons, SCIP_SOL};
 use std::collections::BTreeMap;
 use std::ffi::{c_int, CStr, CString};
 use std::mem::MaybeUninit;
@@ -152,7 +152,7 @@ impl ScipPtr {
         vars
     }
 
-    pub(crate) fn conss(&self) -> Vec<Rc<Constraint>> {
+    pub(crate) fn conss(&self) -> Vec<*mut SCIP_Cons> {
         // NOTE: this method should only be called once per SCIP instance
         let n_conss = self.n_conss();
         let mut conss = Vec::with_capacity(n_conss);
@@ -162,8 +162,7 @@ impl ScipPtr {
             unsafe {
                 ffi::SCIPcaptureCons(self.raw, scip_cons);
             }
-            let cons = Rc::new(Constraint { raw: scip_cons });
-            conss.push(cons);
+            conss.push(scip_cons);
         }
         conss
     }
@@ -279,7 +278,7 @@ impl ScipPtr {
         lhs: f64,
         rhs: f64,
         name: &str,
-    ) -> Result<Constraint, Retcode> {
+    ) -> Result<*mut SCIP_Cons, Retcode> {
         assert_eq!(vars.len(), coefs.len());
         let c_name = CString::new(name).unwrap();
         let mut scip_cons = MaybeUninit::uninit();
@@ -298,7 +297,7 @@ impl ScipPtr {
             scip_call! { ffi::SCIPaddCoefLinear(self.raw, scip_cons, var.raw, coefs[i]) };
         }
         scip_call! { ffi::SCIPaddCons(self.raw, scip_cons) };
-        Ok(Constraint { raw: scip_cons })
+        Ok(scip_cons)
     }
 
     /// Create set partitioning constraint
@@ -306,7 +305,7 @@ impl ScipPtr {
         &self,
         vars: Vec<Rc<Variable>>,
         name: &str,
-    ) -> Result<Constraint, Retcode> {
+    ) -> Result<*mut SCIP_Cons, Retcode> {
         let c_name = CString::new(name).unwrap();
         let mut scip_cons = MaybeUninit::uninit();
         scip_call! { ffi::SCIPcreateConsBasicSetpart(
@@ -321,7 +320,7 @@ impl ScipPtr {
             scip_call! { ffi::SCIPaddCoefSetppc(self.raw, scip_cons, var.raw) };
         }
         scip_call! { ffi::SCIPaddCons(self.raw, scip_cons) };
-        Ok(Constraint { raw: scip_cons })
+        Ok(scip_cons)
     }
 
     /// Create set cover constraint
@@ -329,7 +328,7 @@ impl ScipPtr {
         &self,
         vars: Vec<Rc<Variable>>,
         name: &str,
-    ) -> Result<Constraint, Retcode> {
+    ) -> Result<*mut SCIP_Cons, Retcode> {
         let c_name = CString::new(name).unwrap();
         let mut scip_cons = MaybeUninit::uninit();
         scip_call! { ffi::SCIPcreateConsBasicSetcover(
@@ -344,7 +343,7 @@ impl ScipPtr {
             scip_call! { ffi::SCIPaddCoefSetppc(self.raw, scip_cons, var.raw) };
         }
         scip_call! { ffi::SCIPaddCons(self.raw, scip_cons) };
-        Ok(Constraint { raw: scip_cons })
+        Ok(scip_cons)
     }
 
     pub(crate) fn create_cons_quadratic(
@@ -357,7 +356,7 @@ impl ScipPtr {
         lhs: f64,
         rhs: f64,
         name: &str,
-    ) -> Result<Constraint, Retcode> {
+    ) -> Result<*mut SCIP_Cons, Retcode> {
         assert_eq!(lin_vars.len(), lin_coefs.len());
         assert!(
             lin_vars.len() <= c_int::MAX as usize,
@@ -398,7 +397,7 @@ impl ScipPtr {
 
         let scip_cons = unsafe { scip_cons.assume_init() };
         scip_call! { ffi::SCIPaddCons(self.raw, scip_cons) };
-        Ok(Constraint { raw: scip_cons })
+        Ok(scip_cons)
     }
 
     /// Create set packing constraint
@@ -406,7 +405,7 @@ impl ScipPtr {
         &self,
         vars: Vec<Rc<Variable>>,
         name: &str,
-    ) -> Result<Constraint, Retcode> {
+    ) -> Result<*mut SCIP_Cons, Retcode> {
         let c_name = CString::new(name).unwrap();
         let mut scip_cons = MaybeUninit::uninit();
         scip_call! { ffi::SCIPcreateConsBasicSetpack(
@@ -421,7 +420,7 @@ impl ScipPtr {
             scip_call! { ffi::SCIPaddCoefSetppc(self.raw, scip_cons, var.raw) };
         }
         scip_call! { ffi::SCIPaddCons(self.raw, scip_cons) };
-        Ok(Constraint { raw: scip_cons })
+        Ok(scip_cons)
     }
 
     /// Create cardinality constraint
@@ -430,7 +429,7 @@ impl ScipPtr {
         vars: Vec<Rc<Variable>>,
         cardinality: usize,
         name: &str,
-    ) -> Result<Constraint, Retcode> {
+    ) -> Result<*mut SCIP_Cons, Retcode> {
         let c_name = CString::new(name).unwrap();
         let mut scip_cons = MaybeUninit::uninit();
         scip_call! { ffi::SCIPcreateConsBasicCardinality(
@@ -449,7 +448,7 @@ impl ScipPtr {
         }
         scip_call! { ffi:: SCIPchgCardvalCardinality(self.raw, scip_cons, cardinality as i32) };
         scip_call! { ffi::SCIPaddCons(self.raw, scip_cons) };
-        Ok(Constraint { raw: scip_cons })
+        Ok(scip_cons)
     }
 
     pub(crate) fn create_cons_indicator(
@@ -459,7 +458,7 @@ impl ScipPtr {
         coefs: &mut [f64],
         rhs: f64,
         name: &str,
-    ) -> Result<Constraint, Retcode> {
+    ) -> Result<*mut SCIP_Cons, Retcode> {
         assert_eq!(vars.len(), coefs.len());
         let c_name = CString::new(name).unwrap();
         let mut scip_cons = MaybeUninit::uninit();
@@ -479,7 +478,7 @@ impl ScipPtr {
 
         let scip_cons = unsafe { scip_cons.assume_init() };
         scip_call! { ffi::SCIPaddCons(self.raw, scip_cons) };
-        Ok(Constraint { raw: scip_cons })
+        Ok(scip_cons)
     }
 
     /// Create solution